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Abstract 



This thesis describes the development of DiFX, the first general-purpose software cor- 
relator for radio interferometry, and its use with the Australian Long Baseline Array to 
complete the largest Very Long Baseline Interferometry (VLBI) pulsar astrometry pro- 
gram undertaken to date in the Southern Hemisphere. This two year astrometry program 
has resulted in the measurement of seven new pulsar parallaxes, which has more than tre- 
bled the number of measured VLBI pulsar parallaxes in the Southern Hemisphere. These 
measurements included a determination of the distance and transverse velocity of PSR 
J0437-4715 with better than 1% accuracy, enabling improved tests of General Relativity; 
the first significant measurement of parallax for the famous double pulsar system PSR 
J0737-3039A/B, which will allow tests of General Relativity in this system to proceed to 
the 0.01% level and also offers insights into its formation and high-energy emission; and 
a factor of four revision to the estimated distance of PSR J0630-2834, which had previ- 
ously appeared to possess extremely unusual x-ray emission characteristics. Additionally, 
the ensemble of refined distance and transverse velocity estimates have enabled a widely 
applicable improvement in knowledge of pulsar luminosities in several wavebands and the 
Galactic electron distribution at southern latitudes. Finally, the DiFX software correlator 
developed to enable this science has been extensively tested and verified against three 
existing hardware correlators, and is now an integral part of the upgraded Long Baseline 
Array Major National Research Facility used by astronomers throughout Australia and 
the world; furthermore, it has been selected to facilitate a major upgrade of the world's 
only full-time VLBI instrument, the Very Long Baseline Array operated by the National 
Radio Astronomy Observatory in the US. 
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1 

INTRODUCTION 



1.1 Thesis motivation 

Radio astronomy, in particular radio interferometry and its high resolution sub-branch 
Very Long Baseline Interferometry (VLBI, discussed in detail in Chapter [3|), is a field in 
which advances in instrumentation - driven in this case by developments in consumer and 
industrial electronics - have enabled rapid, ongoing advances in science, by expanding the 
parameter space scientists can explore. This has led to a dependency between engineer and 
scientist which is rarely seen in other fields of astronomy - most radio astronomers have 
at least a passing knowledge of the systems they use, and many are themselves developers 
as well as users of instruments. This is especially true for part-time VLBI arrays such as 
the Australian Long Baseline Array (LBA). 

Another field of study with a strong overlap between engineer and astronomer is pulsar 
astronomy. Radio pulsars (discussed in detail in Chapter [2]) are rapidly rotating neutron 
stars that emit radiation from their magnetic poles. Neutron stars form from the collapsed 
cores of once-massive stars following a supernova explosion. Due to the pulsar's very high 
moment of inertia, the pulsar spin period P is typically very stable. The misalignment of 
the rotation and magnetic axes leads to the radiation being observed as a series of pulses 
(dispersed in frequency by intervening ionised matter) at Earth. Analysis of pulsar data 
typically requires dedicated, high speed signal processing, which has led to most pulsar 
groups developing and deploying their own digital electronic systems on a telescope by 
telescope basis. 

VLBI is an integral tool for the study of pulsars, allowing the determination of kine- 
matic parameters of individual pulsars in a (relatively) precise and model-independent 
fashion. This use of high-resolution observations to accurately measure object positions is 
known as astrometry. As discussed in Section [2321 the addition of independent kinematic 
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Table 1.1. Comparison of major VLBI arrays 



Array 
name 


Array 
stations 


Maximum station 
data rate (Mbps) 


Maximum baseline 
length (km) 


Aetive observing 
(weeks/year) 


VLBA^ 


10 


512 


8600 


52 


EVN 2 


18 


1024 


10000 


10-15 


LBA (S2 - pre-2005) ^ 


6 


128 


1700 


3-4 


LBA (DiFX - post-2005) ^ 


6 


1024 


1700 


3-4 



/ /www. vlba.nrao.edu/ 

/ /www. evlbi.org/intro/ intro.html 

/ /www. atnf.csiro.au /vlbi/ 



IlLLp 

^http: 
^htto: 



information allows the calculation of geometrical effects which alter the observed arrival 
time of pulses. If the signature of the annual orbital parallax on pulsar position imposed 
by the Earth's motion around the Sun can be detected, the resultant determination of pul- 
sar distance can be used to accurately calibrate the pulsar luminosity at all wavebands, 
as well as further refining the pulse arrival time corrections. 

While Southern Hemisphere instrumentation has played a crucial role in the study of 
pulsars - the Parkes and Molonglo radiotelescopes in Australia have discovered over half of 
the known radio pulsar population - few pulsar VLBI observations have b een made from 



the Souther r i Hemisphere. Th r ee pr evious Southern Hemisphere surveys (jPodson et al 



20031; iLeggd, 



2002 



Bailes et al 



1990) have resulted in the measurement of two pulsar par- 



allaxes, whereas 16 Northern Hemispher e parallaxes were pu blished at the time of writing, 
with nine obtained in a single program (jBrisken et al.l . |2002| ) . This is primarily due to the 
capabilities of the American Very Long Baseline Array (VLBA) and the European VLBI 
Network (EVN) instruments, both of which possessed advantages in recording bandwidth, 
support and observation cadence compared to the LBA, which is the only Southern Hemi- 
sphere VLBI array. The LBA, VLBA and EVN antennas are shown in Figure 11.11 - full 
details of these arrays are shown in Table II. 1[ 

Despite the advantages of Northern Hemisphere arrays, there are many unique pulsars 
which lie too far south to be effectively observed from the Northern Hemisphere such as the 
unique double pulsar system PSR J0737-3039A/B, the longest period radio pulsar PSR 
J2144-3933, and the nearest and brightest millisecond pulsar PSR J0437~4715. These 
objects and others offer insights into pulsar formation, evolution and many other related 
fields of research, but have yet to be the subjects of detailed study at the highest angular 
resolution. The LBA offers the only means to undertake high angular resolution studies 
of these objects. 



1.1. Thesis motivation 
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The key impediment faced by the LB A at the outset of this thesis (in early 2005), 
compared to the VLBA and EVN, was a lack of sensitivity. As discussed in Chapter [3l 
the sensitivity of an interferometer is proportional to the square root of the bandwidth of 
the signal it accepts, which is limited by the digital sampling, recording and processing 
hardware employed by the array. As shown in Table [LTl in 2005 the LEA was significantly 
limited in the bandwidth it could record, compared to the VLBA and EVN. This meant 
that targeting the most scientifically desirable Southern Hemisphere pulsars, many of 
which are faint radio sources, would require impossibly large amounts of telescope time to 
obtain sufficient sensitivity. An upgrade of the LEA was thus the only feasible alternative 
to obtain astrometric information on these objects. 

This upgrade involved replacing the existing tape-based recorders and signal processing 
hardware (the " correlator" , discussed in Chapters E] and S]) with disk-based recorders and 
an alternate correlator capable of handling higher data rates. At the commencement 
of this thesis, new disk-based recorders were being tested with a preliminary correlator 



based on a software algorithm running on a small supercomputer (IWestl . l2004l ) . Despite 
verifying the functionality of the disk-based system, this initial "software correlator" was 
too slow (taking weeks to correlate a day's observing) for production usage. A refined, 
more efficient software correlator was required, which was developed during this thesis 
and became known as DiFX (short for Distributed FX correlator - the FX terminology is 
explained in Chapter H]). 

Thus, from the outset, this thesis aimed to address the twin goals of a developing a 
flexible, powerful, and efficient software correlator to make use of the higher bandwidths 
available with a disk-based system, and the integration of that correlator into the LEA 
to produce an instrument with the flexibility and sensitivity necessary to successfully 
undertake an astrometric program encompassing the most scientiflcally fruitful Southern 
Hemisphere pulsars. The undertaking of this astrometric program also necessitated the 
development of signiflcant new algorithms and tools for astrometric data reduction, and 
the characterisation and improvement of many areas of LEA operations. 

While the development of the DiFX software correlator was a necessary precursor to 
the desired astrometric science, the improved sensitivity and flexibility of the new system 
extended the capabilities of the LEA for all science targets. Indeed, DiFX has also been 
adopted by several new or upgraded arrays external to Australia, most notably the VLEA. 
Section 14.61 briefly discusses some science highlights external to this thesis that have been 
obtained using the DiFX software correlator on both the LEA and other arrays. 
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Figure 1.1 The location of antennas regularly participating in VLB A (red), EVN (blue) 
and LB A (green) observations. Antennas which are sometimes added to one or more arrays 
on an ad-hoc basis, or that belong to other arrays such as the Japanese VLBI Network 
(JVN) or the Korean VLBI Network (KVN) are shown in black. 



1.2. Thesis outline 
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1.2 Thesis outline 

An overview of the historical studies and present understanding of pulsars is given in 
Chapter [2l along with a discussion of the pulsar science which is possible through the 
use of VLBI astrometry. Chapter [3] presents a conceptual and mathematical overview of 
radio interferometry, including VLBI, and covers the application of VLBI to astrometric 
observations. Chapter H] covers the development, testing and verification of DiFX, the 
final version of the software correlator developed to fulfil the first primary goal of this the- 
sis. Chapter [5] examines the post-correlation data analysis undertaken on all astrometric 
datasets, showing the transformation from correlated data to pulsar positions at a given 
epoch. Chapter [6] highlights the results obtained from the LBA astrometric program and 
shows the implications of the measured pulsar distances and kinematics, both for each 
pulsar individually and for population studies as a whole. Concluding remarks are made 
in Chapter [71 



2 

PULSARS 



2.1 Discovery and studies 



When the existence of highly compressed stellar objects consisting primarily of neu- 
trons - neutron st ars - wa s postulated as a possible result of a supernova explosion by 
Baade and Zwickvl in Il934l . the field of radio astronomy was barely taking its first tenta- 
tive steps. A third of a century later, however, it would be radio astronomy that provided 
a remarkable confirmation of the existence of neutron stars, beginning with PhD student 
Jocelyn Bell noticing a periodic "little bit of scrufF' while obse rving at a Canibridg e ra- 
diotelescope. This discovery was published the following year (iHewish et al.l . Il968l ) and 
short l y afte r the rotating neut ron star origin of the signal was independently proposed by 



Gold 



andlPacinil (|l968l l. 



These periodic radio sources became known as "pulsars", and a flood of theoretical and 
observational results followed the initia l disc overy. Before the end of the dec ade, pulsars 

19691 ) and op tical (ICocke et al.l . I1969I 1 wavebands. 



were detected in the x-ray (iFritz et al 



with detection in gamma rays following shortly after (jFazio et al.l . ll972l ). By 1980 over 300 
pulsars had been discovered and the neutron star origin, with the radio emission powered 
by the conversion of rotational energy, was well established. 

Since that time, however, a series of observational surprises have shown that the neu- 
tron stars can manifest themselves in a variety of guises: 



recycled or millisecond pulsars, with lower (~ 10*^ gauss) magnetic field strengths 
and spin freque ncies in the hundr e ds of Hz, formed throu gh accretion of matter in a 
binary system (lAlpar et al.l . Il982l : 



van den Heuvel 



197,^ ): 
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Anomalous X-ray Pulsars ( AXPs) , which emit more energy than can be explained by 
their spindown rate, and Soft Gamma Repeaters (SGRs), both of which are believed 
to be magnetars - neutron stars with extremely high (> 10^^ gauss) magnetic fields, 
where the m agnetic field decay powers repe ated powerful outbursts in x-rays and 
gamma-rays (jThompson and Duncanl . Il996l ) ; and 



Nulling pulsars, intermittent pulsars and Rotating Radio Transients (RRAT s ), where 



the radio emission is intermitt ently suppressed (see e.g. IWang et al 



197C 



2003; 



Backed . 



McLaughlin et al.l . |2006| ). In the case of RRATs, as few as one pulse in thou- 



sands is emitted. 



As the objects studied in this thesis are all rotation-powered, non-nulling radio pulsars, 
the remainder of this chapter will focus on these objects. 



2.2 Current understanding 



2.2.1 Formation 

Ordinary radio pulsars are believed to form in the supernova explosions which result when 
massive stars exhaust their supply of the light elements which had fueled nuclear fusion. 
With the abrupt removal of radiation pressure which had supported the star against grav- 
ity, a rapid contraction follows. Depending on the stellar core mass, one of three compact 
objects is formed in the final contraction - a white dwarf, neutron star or black hole. With 
a core mass of less than roughly 1.4 solar masses (Mq), the stellar material becomes com- 
pletely ionised during the core collapse, and the Fermi pressure of the resultant degenerate 
electron gas grows until, when the core is several thousand km in radius, it balances the 
gravitational force and a stable, cooling white dwarf remains. However, for cores exceed- 
ing the Chandrasekhar limit of roughly 1.4 Mq, the rising degenerate electron pressure is 
not sufficient to halt further collapse into a denser state - a process first recognised by 



Chandrasekhad (Il93ll ). This violent contraction and ejection of stellar material is known 
as a core collapse supernova. 

For core masses exceeding the Chandrasekhar limit, the collapsing material reaches 
densities and temperatures sufficient to fuse electrons and protons into neutrons and elec- 
tron neutrinos via the process of inverse beta decay: 

e~+p^^n + Ve (2.1) 



2.2. Current understanding 
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The escape of these neutrinos from the collapsing core cools the collapsing material, 
and the simultaneous loss of thermal and Fermi pressure removes any means for the stellar 
material to resist gravitational collapse. It should be noted that the timescale of the neu- 
trino emission remains somewhat uncertain, due to the difficultv of s i mulat ing the extreme 



2002. and references 



environment of the supernova collapse (see e.g. iFrver and Warrenl . 
therein). The collapsing matter largely conserves angular momentum and magnetic flux, 
and thus the initially modest rotation speeds and magnetic fields of the progenitor star 
are amplified immensely as the core compresses. 

The production of a neutron star or black hole depends upon the core mass - for 
masses below about 3 Mq, the Fermi pressure of the degenerate neutron fluid grows until 
it balances the gravitational pressure at a stellar radius of roughly 10 km. The resultant 
neutron star has a core density of > 10^^ g cm^, several times denser than an atomic 
nucleus. The inferred composition of neutron stars is discussed further in Section r2.2.2i For 
collapsing core masses above about 3 Mq , even neutron degeneracy pressure is insufficient 
to stop the collapse, and the core is pred i cted t o collapse completely to form a black hole 
in a hypernova explosion ( Iwamoto et al.l . 1998 ). 

Observations of pulsars have shown that the simple core collapse model described 
above alone cannot completely explain typical pulsar characteristics. One of the chief 
problems is the extremely high space ve l ocity which many pulsars have been observed to 
possess. Recent estimates (jHobbs et al.l . l2005l ) put that average pulsar 3D birth velocity 
at 400 km s~^, with the fastest known pulsar (PSR 61508+55; IChatterjee et al.l . |2005| ) 
possessing an astonishing transverse velocity of 1100 km s^^! These velocity values are 
much higher than thos e possessed by the massive sta rs which are neutron star progenitors 



(~ 20 km s ^; see e.g. lFeast and Shuttleworthl . Il965l ). which along with the small number 



of pulsars in binary systems implies that some physica l process imparts a large velocity 
on most neutron stars at birth (e.g. iDewev and Corded . 119871 : iBailed . 1 19891 ) . Whilst the 
disruption of binary systems may account for some pulsar velocities, it appears that some 
kind of '"kick" mechanism during the formation process is required to adequately explain 
the full range of observed systems. Counter-examples, such as the PSR J0737-3039A/B 
system discussed in Section I6.1.2|, seem to imply that kicks are not universal, further 
complicating interpretations of the physical mechanism. 



10 



Chapter 2. PULSARS 



While many theories have been advanced, generally requiring an asymmetry in the 
collapse and / or neutrino emission during the su pernova formation, or asymm etric elec- 
tromagnetic radiation after the collapse (see e.g. iFrveii . l2004l : iLai et al.l . |200l|), the exact 
nature of the kick mechanism remains unclear. It is important to note that many theories 
of pulsar kicks predict that the kick is aligned with the pulsar spin axis, w hich is tested ob 
serva t ionally by mea suring the polarisation position angle of pulsars (e.g. 



2005; 


Rankin 




2007) 


(e.g. 


Gaensler et al. 



Johnston et al 



20071) or th e position angle of an observed pulsar wind nebula and/or jet 

and comparing to the velocity position 



20021 : 



Helfand et al. 



2001 



angle. Thus, studies of the space velocities of pulsars allow important insights into their 
formation processes. 



2.2.2 Composition 

The composition of the end state of matter in the interior of a neutron star is still the 
subject of controversy. The Equation of State (EoS) of the material at the core of a neutron 
star, which describes the relationship between density and pressure, is a much sought-after 
result which could be obtained from a simultaneous measurement of a neutron star mass 
and radius (see eg. 



Lattimer and Prakash 



20071 ). The extreme environn ient mea,ns tha t 



exotic forms of matter could exist in the core, such as unconfined quarks (jPrakashl . 120071 ). 
The difficulty of such measurements means that to date, a wide range of EoS's remain 
permitted. The discussion below assumes neutron stars do not contain exotic material. 

Figure 12.11 shows the generally accepted taxonomy of a "normal" neutron star. It can 
be broadly divided into the "crust" and "core" regions - although the phase transitions 
between the regions are poorly understood and intermediate layers could exist. The crust 
consists of atomic nuclei and free electrons, since the pressure near the neutron star surface 
is low enough to permit nuclei to remain intact. The atomic nuclei are locked in a solid 
lattice-like structure. This rigid crust "freezes" the magnetic field configuration of the 
neutron star in place. As the density increases further from the surface, free neutrons 
(whi ch are predicted to exh ibit superfiuidity) are present along with nuclei and electrons 
(e.g. 



Sandulescu et al 



2004 ) - this is shown as the "inner crust" in Figure 12. 1[ At a depth 
of several km, the pressure becomes too great for atomic nuclei to exist and a transition 
to a neutron-only environment occurs. 
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Figure 2.1 The interior structure of a neutron star. The (predominantly iron) nuclei in 
the crust form a solid lattice, "freezing" the magnetic field of the neutron star in place, 
while free neutrons in the inner crust and core are believed to exhibit superfluidity. 



The neutron core of a neutron star is expected to be superfluid, although the math- 
ematical treatment of the formation of Cooper pairs of neutrons at nuclear density is 
exceedingly complex. Some free protons are also expected to exist in the core, forming a 
superconducting dynamo which is the source of the neutron star's magnetic field. Obser- 
vational support for this model has come from pulsar glitches - events where the rotational 
frequency undergoes a sudden increase (spin-up), before recovering over a longer period 
to resume the steady spin-down caused by the loss of rotational energy. The spin-up 
can be explained by a transfer of angular momentum between the crust and core facil- 
i tated by vortices in th e superfluid core, changing the neutron star's moment of inertia 
([Larson and Linkl . 12002 ). 
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2.2.3 Emission 



Despite decades of intense study, a complete picture of the pulsar emission mechanism has 
proved elusive. However, some aspects are well understood, and the generally agreed facts 
are presented below. 

A schematic diagram of an ordinary radio pulsar is shown in Figure 12.21 The "light 
cylinder" reflects the distance from the pulsar at which a particle with the same angular 
velocity as the pulsar would be required to travel at the speed of light - thus, the region 
outside the light cylinder is prevented from co-rotating with the pulsar. As shown in 
Figure 12.11 the solid crust of the neutron star locks the magnetic field of the pulsar and 
forces it to co-rotate with the star. 

The pulsar magnetosphere (shown in Figure 12. 2p is the relativistic ch a rged plasma 
which co-rotates with the pulsar, first postulated bv iGoldreich and JulianI (|l969l ). This 
co-rotation of the magnetosphere means that magnetic field lines which originate close to 
the magnetic axis of the pulsar are forced to remain open, since their closure would require 
them to cross the light cylinder, meaning the plasma locked to these field lines would be 
traveling faster than the speed of light. Thus, for some region around each magnetic pole 
the field lines cannot close - these regions are known as the "polar caps" . 

The presence of charged plasma in the pulsar magnetosphere can explain both the 
non-thermal and thermal emission of pulsars - the non-thermal emission by synchrotron 
emission from electrons and positrons spiralling away from the pulsar around the open 
magnetic field lines, and the thermal emission from the polar cap regions, which would be 
heated by the impact of infalling relativistic material. The charged particles which form 
the pulsar magnetospher e are c ontinually replenished through a pair-production process 
(jPaughertv and Hardina . Il982l ). fed by high-energy 7-ray photons, themselves produced 
by curvature radiation from particles accelerated along the curved magnetic field lines in a 
very large electric potential produced in a vacuum gap somewhere in the magnetosphere. 
However , the process which p roduces the coherent radio emission is still not understood 
(see e.g. 



Lyutikov et al. 



I999I . and references therein). 



The precise site of the massive electric potential required t o produce very high energy 
ray photons is not yet well understood. Early models such as 



Ruderman and Sutherland 



(|l975l ) proposed the site of th e acceleratiori was deep in the magnetosphere, ii ear the polar 



cap, while later models (e.g. ICheng et al 



198i 



Chiang and Romanil . Il994l ) proposed a 



location much higher in the magnetosphere, in a region known as the outer gap. The 
difficulty of testing these models with the available observations means that one or both 
could be correct, and the location of the gap could vary from pulsar to pulsar. 
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Rotation Axis 




Figure 2.2 Representation of the pulsar magnetosphere. Magnetic field lines with an open- 
ing angle greater than a critical value are forced to remain open by the co-rot ating charged 



plasma locked to closed field lines within the light cyclinder (jMichell . Il974l ). The resul- 
tant rotating magnetic dipole emits electromagnetic radiation, and pair production in the 
intense magnetic field provides a source of charged particles which are either accelerated 
away from the star - the "pulsar wind" - or back onto the polar cap (the region defined 
by open field lines), heating it. The heated polar cap is believed to be a source of thermal 
x-rays. 
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Another avenue for understanding the pulsar emission mechanism is high er-frequency 



obser vations. Rotation -powered pulsars 



l ave be en detected in the optical (e.g. 



20041 ) ■ ultraviolet (e.g. iKargaltsev et al 



Zharikov et al. 



20041 ). and x-ray (e.g. IKargaltsev et al.l . l2006l ) 



wavebands, frequently with contradictory results. While some observations have sug- 
gested a non-thermal spectrum for higher energy radiation, presumed to originate in the 
pulsar magnetosphere, others have found emission consistent with a solely thermal source. 
Thermal emission will be generated over the entire neutron star surface, but as shown 
in Figure 12.21 the polar cap regions are expected to be heated relative to the remainder 
of the surface and may dominate the high-energy thermal emission. Models of pulsar 
surface temperatures are presently poorly constrained, and hence accurate luminosities 
for rotation-powered pulsars detected in the optical to x-ray wavebands are essential for 
ongoing attempts to understand the pulsar emission mechanism. The chief uncertainty in 
many pulsar luminosity estimates is the large uncertainty in the pulsar distance. Remov- 
ing these uncertainties for specific objects through direct measurement of pulsar distances 
is one of the applications of this thesis. 

As shown in Figure 12.21 the misalignment of the rotational and magnetic axis of the 
pulsar leads to the emission beam tracing a conical shape on the sky. If the beam of 
radiation intersects with Earth during its path, an observer on Earth can detect periodic 
pulsed emission. 



2.2.4 Isolated pulsar evolution 

For ordinary radio pulsars, the emission of electromagnetic radiation due to the rotating 
magnetic dipole, governed by the classical electrodynamics, is the primary mechanism for 
energy loss. This "magnetic braking" leads to a steady increase in the pulsar's rotational 
period as rotational energy is lost. Under the assumption that other forms of energy loss 
are negligible, the followi ng relation between P, P a nd surface magnetic field strength B 



can be obtained (see e.g. iLorimer and Krameij . l2005l ): 



B = \l^ J ^ PP (2.2) 
V 87r2 i?6 sin2 a ^ ' 

where R is the neutron star radius (usually estimated as 10 km), I is its moment of inertia 
(usually estimated as 10^^ g cm^) and a is the angle between the magnetic moment of the 
neutron star and its spin axis. 

A common way of visualising the known pulsar population is to plot logP against 
log P, or equivalently (under the assumption of pure magnetic dipole braking) log P against 
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log B using equation 12.21 This is shown in Figure 12.31 This diagram shows three distinct 
pulsar populations - "normal" radio pulsars in the centre, with 10^^ < B < 10^^ G and 
10^^ < P < 10^ s, recycled pulsars with B < W^^ G and short rotational periods, and 
magnetars/AXPs with B > 10^^ G and generally long rotational periods. 

While the evolution of the magnetic field strength of non-recycled, rotat ion-powered 



pulsars over their lif e times has been the subject of consi derable debate (see e . g.lGoldreich and Reiseneggeii . 



19921 : iBhattacharval . |2002| ) . numerical simulations (e.g. iBhattacharva et al.l . Il992l ) support 



the view that the field strength does not evolve with time, and so constant field strength 
is generally assumed. Additionally, pulsars are often assumed to be born with initial spin 
periods of order 1-20 milliseconds, which is reasonable based upoi i conservat i on of angular 
momentum of the progenitor, and is backed up by simulations (jOtt et al.l . 120061 ). How- 
ever, various observational results suggest that pulsa rs can also be bor n with considerably 
longer spin periods (e.g. P SR J1811-1925, 65 ms, borii et al.l \l99^ : PSR J0538+2817, 
2OO3I ) . Nevertheless, generally speaking, pulsars are born close to 



140 ms, 



Kramer et al. 



the left of the P-P diagram and evolve towards the right along lines of constant B. The 
pulsar "characteristic age" Tc can be estimated by assuming the initial spin period is 
negligible compared to the current period using: 



P 



(n- 1)P 



(2.3) 



where n is the braking index, which is equal to 3 for pure magnetic dipole braking in a 
vacuum. 

The pulsar's true age could vary considerably from its characteristic age if the initial 
spin period Pq was a significant fraction of the current spin period, or if the assumption 
that magnetic braking is the dominant form of energy loss is incorrect. An alternative 
way of estimating pulsar ages is to use the "kinematic" age, which is calculated when the 
pulsar's position, proper motion and birth location are known. This method can generally 
only be used for pulsars associated with a known supernova remnant, or for whom the 
uncertainty in birth location (pulsars are generally assumed to be born close to the Galac- 



Cordes and Chernofi 



19981 . 



tic plane, with a scale height of approximately 100-130 pc 

Faucher-Giguere and Kaspil2006 ) is unimportant, such as pulsars with a Galactic height 



many times the scale height, with a well-determined vertical proper motion. Comparison 
of kinematic and characteristic ages can lead to constraints on the birth location or initial 
spin periods of pulsars. 

It is apparent from Figure 12.31 that a region of parameter space with low magnetic 
field strength and long periods is completely depopulated of pulsars. This pulsar "death 
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Figure 2.3 Plot of known pulsars showing the distribution of spin period P and surface 
magnetic field strength B (as calculated from P using Equation I2.2p . Binary systems are 
shown as semi-filled circles and isolated neutron stars with dots. Magnetars lie in the 
upper right-hand corner of the plot, and recycled pulsars in the lower left-hand corner. 



zone" arises naturally from the pair-production cascade emission mechanism discussed in 
Section r2.2.3l - once the curvature radiation no longer produces photons with sufficient en- 
ergy to initiate a pair-production cascade, the observable radio emission ceases. However, 
several observed pulsars (including most notably the P = 8.5 s PSR J2144-3933) contra- 
dict this simple picture, being observable despite being past the pulsar death line. The 
existence of observable long-period pulsars such as PSR J2144-3933 may be explained by 
invoking inverse Compto n scattering as the source of the pair-production cascade, replac- 
ing curvature radiation (jZhang et al.l . I2OOOI ). Observations of pulsars in the death zone 
which allow some insight into their fundamental properties such as luminosity may help 
to resolve the questions as to how their radio emission is maintained. 



2.2.5 Binary pulsars 

Pulsars exist in a wide range of binary systems, orbiting other neutron stars (including 
the famous double pulsar system PSR J0737-3039A/B), white dwarfs, main-sequence 
and post main-sequence stars. Their presence in binary systems is a prerequisite for the 
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formation of recycled pulsars, where the accretion of matter from a companion star that 
is overflowing its Roche lob e transfers angular momentum to the pul sar, spinning it up to 
millisecond periods (see e.g. iBhattacharva and van den Heuvell . Il99ll . for a comprehensive 
discussion). Although rarer than lone pulsars, binary pulsars offer a number of unique 
science opportunities, discussed below. 

The kinematics of binary systems including a pulsar offer an insight into the supernova 
events which form neutron stars, since the binary system must survive the supernova 
undisruptecill. The orbital semi-major axes and eccentricities of known binary pulsars 
allows some constraints to be placed on the types of progenitor stars and supernovae that 
lead to these systems. 

Pulsars in close binary orbits travel at relativistic speeds and offer the possibility to 
test the predictions of General Relativity (GR) against alternate theories of gravity in the 
strong-field gravity regime. Examples of post-Keplerian effect include decay in orbital 
period Ph due to gra vitational wave emission, relativistic orbital precession and Shapiro 
delay (IShapird . Il964l ) . The detection of these effects is dependent on the precision timing 
of pulse arrival, which is discussed in Section [2.3.11 

Finally, mergers of compact objects in binary systems are e xpected to be one of th e 
first sources of gravitational waves to be directly detected (e.g. iBelczynski et al.l . 12002). 
The population size of relativistic binaries in our Galaxy is crucial when estimating the 
frequency of merger events throughout the local Universe, and thus the probability of 
success for the Laser Interferometer Gravitational Wave Observatory (LIGO). Estimates 
of population size for relativistic binaries depend on the spatial density and luminosity 
function of the systems, which require accurate distances. Thus, observations of existing 
binary systems can contribute to the expected frequency of gravitational wave events, and 
hence event detection rates with LIGO. 



2.3 Observing pulsars in the radio waveband 



As shown in Section [2.2.31 pulsars generate beams of coherent, broad-band radio emission 
which is observed as a pulse train at Earth, due to the pulsar's rotation sweeping the beam 
past Earth. Pulsars are generally observed to have st eep spectra - the mean spectral index 
for normal radio pulsars is —1.8 (jMaron et al.l . l200d ). 

In order for the pulsar signal to propagate to Earth, it must pass through the pulsar's 



^Globular cluster binaries, which have a much higher incidence of interactions, do not necessarily offer 
the same insight 

^ Those not predicted by classical mechanics 
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Figure 2.4 The effect on pulsar radiation of travelling through a medium of non-zero 
density. The broadband pulses, which are initially aligned in frequency, are dispersed by 
ionised material along the line of sight. Density variations cause refractive and diffr active 
scintillation and scattering. 



local environment, the interstellar medium (ISM) and Earth's own atmosphere. Each of 
these environments is typically composed, at least in part, by non-uniformly distributed 
ionised matter which interacts with the radio waves. Essentially, the radiation traverses a 
path of continually varying refractive index, which causes dispersion, scintillation (both re- 
fractive and diffractive) , scattering and (for polarized radiation in the presence of magnetic 
fields) Faraday rotation. These effects are summarised below in Figure [231 



Much of the unique science made possible by pulsars depends upon their intrinsic 
rotational stability, enabling their pulsed signals to be taken as accurate clock ticks. For 
this approach to be viable, the propagation effects discussed above must be overcome, 
along with a host of other error sources. This discipline of pulsar timing is discussed 
below. 
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2.3.1 Pulsar timing 



Pulsar timing determines a pulse time of arrival (TOA) by cross-correlating the observed 
pulse profile from an observation (obtained by averaging in time and frequency) and a 
template. This is then compared with a timing model, and a "timing residual" obtained. 
A bootstrap procedure follows, with the model undergoing refinement until an optimal 
model is obtained. Typical timing residuals can be very small fractions of a pulse period 
(e.g. 0.003% for PSR J0437-4715; Iverbiest et al.l . boosll. For most pulsars, the intrinsic 



Hot an et al 



average pulse profile is very stable over time (see e.g. 

is not universal. Geodetic precision can cause secular change s in pulse profi^ 



2004 ). although this 
e over long 



timescales, an effe ct which has been seen in PSR B1913+16 (jWeisberg et al, 



PSR J1141-6545 (Hotanetal 



19891 ) and 



20051 ). On i nuch shorter timescales, so -called "mode- 



Liu et al. 


2006 


Lvne, 


1971 



changing" pulsars such as PSR B0329+54 (e.g. lLiu et al.l . l2006l : lLvnd . ll97ll ) switch between 
two or more modes, in which pulse components vary in relative and absolute strength. Such 
complications are not relevent to this thesis and are not considered further. 

The averaging of recorded data in time requires an accurate model for pulse arrival 
times - the pulsar ephemeris. The instantaneous position of the source (pulsar) and 
observer (telescope on Earth) must be calculated to a high degree of precision. The pulsar 
reference position, proper motion, binary motion (if applicable), rotational period and 
period derivative must be known, as well as the Earth's ephemeris and telescope location. 
This a priori model is summarised below in Figure [231 

Averaging the received pulsar signal in frequency requires the removal of the dispersive 
effects of the ISM. An example of a dispersed pulsar signal is shown in Figure 12.61 The 



time delay in seconds experienced by a pulse at frequency z/ GHz can be expressed as: 

DM 



Td 



2.41 X 102zy2 



(2.4) 



where DM is the so-called "Dispersion Measure" associated with a pulsar. DM is defined 
as the integral of electron column density along the line of sight to the pulsar, quoted here 
in pc cm~^. Observed values of DM range from < 5 pc cm^'^ for very nearby pulsars, to 
> 1000 pc cm~'^ for distant pulsars in the Galactic plane. While DM is often assumed to 
be constant, the relative motion of the pulsar. Earth, and ISM actually lead to continual 
small changes in DM due to the changing electron content along the line of sight, and for 
precision timing the time variation of DM must be measured and applied. 
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Proper motion 




telescope ephemeris 



Figure 2.5 Components of the geometric model used in recording pulse arrival times. The 
orbital and transverse velocity and acceleration of both the pulsar and observer must be 
modeled, requiring an ephemeris for both the pulsar and Earth. 
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9 9.4 M 



Figure 2.6 The intensity of pulsar PSR J0437~4715, shown as a function of frequency and 
pulse phase (time modulo the pulse period, expressed in units of fractional pulse period). 
The increasing delay of the signal with decreasing frequency is clearly apparent. 



For pulsar timing, this frequency-dependent dispersion is simply an inconvenience to 
be characterised and removed, as discussed below. However, it also provides an accurate 
measure of the ionised content of the ISM lying between the pulsar and observer, which 
can be translated into an estimate of the pulsar distance, given an estimation of the density 
of the ionised ISM along the line of sight. Widely used models of the Galactic electron 
distribution, which allow calculation of the io n ised I SM content along arbitrary sightlines, 
have been constru cted by Tavlor and Cordes ( IQoj ). which is hereafter referred to as the 
TC93 model, and 



Cordes and Lazio 



liereafter referred to as the NE2001 model. 



Since pulsar luminosities vary over many orders of magnitude, this provides the most useful 
estimate of distance which can be obtained for the entire pulsar population. However, since 
the density of the ionised ISM can also vary over many orders of magnitudes on small scales, 
feedback into the electron distribution models in the form of model-independent distances 
is crucial to improve the quality of distance estimates for the bulk of known pulsars. 
Methods of obtaining such model-independent distances through VLBI are discussed in 
Section f2.3.2l below, and demonstrated in Chapter El 

Removal of the frequency-dependent delay from the observed pulsar signal, to allow 

to noise ration (SNR), can be 



the summation of data in frequency to improve the signa 
accom plished in one of two ways. Incoherent dedispersion (jVoute et al 



2002 



Large et al 



19691 ) makes use of a filterbank to divide the observed radio band into narrow frequency 



channels, and compensates the delay on a channel by channel basis, with delays appropriate 
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for the mean channel frequency. Since the channels remain a finite width MHz, there 
is some small residual smearing which can be calculated for an observing frequency of v 
GHz by: 



8.3 X X DM 



/iS 



(2.5) 



Alternatively, coherent dedispersion may be employed. Essentially, this approach ap- 
plies a suitable filter to the baseband data (containing the inverse of the transfer function 
of the ISM) before the channelisation process, minimising finite bandwidth effects (how- 



ever, the fundamenta 
was first suggested by 
in high-precision timing campaigns (eg 



limitation of the original sampling time remains). This approach 
Hankins and RickettI (119751). a.nd is becoming increa.singly prevalent 



Hotan et al. 



20061; 



Hessels et al 



20061). 



As shown in Figure 12. 4^ refractive scintillation is caused by large-scale structure in 
the ISM, which acts as a large lens focussing or defocussing the pulsar radiation. This 
naturally leads to amplitude fluctuations in the pulsar signal, and since it is caused by 
large-scale structure, acts over long time periods and large observing bandwidths. Diffrac- 
tive scintillation, on the other hand, is caused by small-scale fluctuations in the ISM, with 
diffraction producing an interference pattern on a plane which the Earth traverses. Vari- 
ations in amplitude are seen as the Earth passes through the diffraction "scintles" due to 
its transverse velocity, and the diffraction pattern itself moves a t the relative speed o f the 



Cordes et al. 



(119861 ) give 



pulsar compared to the ISM where the diffraction is occurring, 
a more detailed overview of the physics of scintillation. 

Since the pulsar velocity is usually much larger than that of the Earth or ISM, scin- 
tillation observations of pulsars can be used to make an estimate of the magnitude of the 
pulsar transverse velocity. This requires an estimation of the pulsar distance and a number 
of simplifying assumptions about the nature of the scattering material. A comprehensive 
discussion of the use of sc intillation studies to make velocity estimates can be found in 
Cordes and RickettI (jl998l ). For the commonly assumed case of a single, do minant thin 



scattering screen, the scintillation speed viss is given (e.g. iGupta et al.l . 11994 ) by: 



viss 



A., 



(2.6) 



where A^ is a constant related to the structure function of the ISM (equ al to 3.85 x 



10^ km s ^ for Kolmogorov tu rbulence in the thin-screen approximation; iGupta et al 



19941 : 



Cordes and Rickett 



19981 ). D is the distance to the pulsar in kpc, X is the ratio of the 



Earth-screen distance to the pulsar-screen distance, is the observing frequency in GHz 
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and Ai'd and At are the decorrelation bandwidth in MHz and decorrelation time in seconds 
respectively. The decorrelation bandwidth and time are determined observationally by 
averaging diffraction scintles and determining the mean bandwidth and time required 
for the pulsar intensity to fall to 1/e of the peak value. The magnitud e of the true 
pulsar velocity \vt\ is related to the scintillation speed by X x \vt\ = Viss (jGupta et al 



19941 : 



Cordes and Rickett 



19981 ). Comparison of predicted pulsar scintillation velocities 



to observed values from VLBI and pulsar timing are made in Chapter [6l 

Inhomogeneities in the ISM also lead to scattering, where reflected "echoes" of the 
pulsed emission are seen after a time delay, as shown in Figure 12.41 Scattering scales 
strongly with frequency, but the exact form of the scaling depends on the distribution of 
material in the intervening ISM - for the commonly assumed Kolmogor ov model of tur- 
bulen ce in the ISM, the frequency dependence of scattering is (e.g. iLee and Jokipiil . 



19751). 



All of these time-variable propagation effects lead to variations in th e pulse arrival tirn e 



estimates. The dominant effect is that of DM variations, as shown by 



You et al 



(120071). 



Scattering variations are less noticeable, and variations due to refractive and diffractive 
scintillation have traditionally been neglected, although simu lations suggest that the ir 
effects are detectable at low frequencies for well-timed pulsars (jFoster and Cordesl . Il990l ). 

A final source of arrival time errors can be instrumental in nature. The propagation of 
signals through analog or digital filterbank and sampling system must invariably involve 
delays, which vary from instrument to instrument, and telescope to telescope. Changes to 
the signal path before the pulsar hardware can also affect instrumental delays. For long 
time series of pulsar observations, which generally span multiple instruments, calibration 
of the unknown relative instrumental delays introduces additional free parameters to the 
timing model. 

A pulsar timing campaign requires regular observations to obtain a series of residual 
delays. Whilst the arrival time errors due to finite signal to noise should be zero-mean, 
Gaussian distributed random noise, incorrectly modeled or neglected effects manifest as 
clear trends in the residual errors. For example, an error in P will result in residuals 
which linearly diverge from 0, while an unmodeled binary companion will lead to periodic 
movement in residuals at the binary period of the system. Historically, analysis of pulsar 
timing data and the fitting of pulsar parameters has used the software package TEMPoR , 



although more recently a more advanced package known as TEMPO^J (iHobbs et al 



is now available, and incorporates support for higher precision experiments than TEMPO, 



^ http : // www.atnf.csiro.au / research /pulsar / tempo / 

^http: / / www.atnf.csiro.au / research /pulsar / psrtime / tempo2 / 
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as well as simultaneous timing of multiple pulsars. 

A large proportion of the exciting science possible using pulsar timing involves bi- 
nary pulsars. For example, exploration of GR effects is generally only possible in binary 
systemjfl. The relevant equations in which a VLBI measurement of kinematics can con- 
tribute to the precision of a timing result are presented below - for an excellent review of all 
the eq uations relevant to pulsar timing, see the pulsar handbook of 



Lorimer and Kramer 



(120051) 



For this thesis, the important timing equations are those dealing with orbital motion 
in binary pulsars. Equation 12.81 shows the factors which contribute to an observed change 
in binary period P^^: 



pobs ^ (2.7) 

pGR ^ pdrag ^ ptid ^ pmA _ /pace ^ pShk\ (2.8) 



where the intrinsic contributions to P^^'^ due to energy loss from the system (P™*) consist 
of relativistic effects such as the emission of gravitation radiation (P^^), atmospheric 
drag (P^^^^), mass loss (-P™^), and tidal dissipation {P^^^'^), and the kinematic contribution 
consists of the relative acceleration of the pulsar to the timing reference point (the 
solar system barycentre) P^^'^ and the apparent acceler a.tion caused by th e pulsar's proper 
motion -Pj^^'^, which is known as the Shklovskii effect ( Shklovskii . 197ol ). The kinematic 
contributions to P^^^ can be expressed as: 



pace 1 



and 



BP ■ (^psr - a bar) (2.9) 

c 



pShk 2 

^ - (2.10) 



Pb cd 



where BP is a unit vector from the solar system barycentre to the pulsar, 7? psr is the 
acceleration of the pulsar system, "a* bar is the acceleration of the solar system barycentre, 
vt is the transverse velocity of the pulsar, d is the distance from the pulsar to Earth 
and c is the speed of light. The acceleration terms "(?bar and 7? psr incorporate Galactic 
rotation, the vertical potential of the Galaxy (and the parent cluster for globular cluster 
pulsars), and any unmodeled nearby perturbing massive bodies. These apparent and 
actual accelerations due to kinematic effects also affect the pulsar's spin period P in a 
similar manner, for both isolated and binary pulsars. 



""A counter-example is microlensing of pulsars, although this is yet to be observed 
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Equations 12.81 - 12.101 show that in the presence of accurate timing data for a binary 
pulsar and sufficiently accurate modehng of one set of parameters, o ne of the contributing 
factor s to Pb can be calculated. This was famously demonstrated bv lTavlor and Weisberg 
( 1983 ) for ^ PSR B1913+16 system. Combining the measured value for Pb 

with the observed and modeled components can yield a limit on the anomalous ac- 
celeration of the pulsar with respect to the Solar System. This anomalous accelera- 
tion can be interpret ed as a change in the value of Newton's gravitational constant G 



( Verbiest et al 



200811. the presence o f a dis tant, massive planetary companion in the So- 



lar System (IZakamska and Tremaind. l2005ll. or an error in the estimate of the Galactic 
gravitational potential ( Bell and Bailesl ." 1996 ). 

Often, the most uncertain contribution to P^^^ is the Shklovskii term i^'^'^, due to the 
large uncertainty which is generally associated with most pulsar distance measurements, 
which contributes a large uncertair ity to the tran s verse velocity. Even in the case of 



excellent timing data, such as the 



Verbiest et al 



(120081) measurement of PSR J0437- 



4715, direct measurement of distance through detection of the annual parallax yields 
an uncertainty of 8%, resulting in an unacceptably large error in ib^'^j if the aim is to 



constrain another contributior 



to P^^^. If the transverse velocity of a pulsar can be 



supplied through an independent measurement of parallax and proper motion, such as 
VLBI astrometry, the uncertainty in P^^^ can be reduced sufficiently to allow significant 
constraints on other terms in P^^*^. An example of this is shown, with limits on G and 
massive planetary companions, for PSR J0437-4715 in Section [6. 1.1[ 



2.3.2 VLBI pulsar observations 

For the purposes of VLBI, pulsars can be considered as unremarkable radio sources, except 
for the possibility of improving the sensitivity of observations through pulsar "gating": 
blanking the telescope data at times when the pulsar flux is low or zero. This has the 
effect of eliminating noise which would otherwise be accumulated during these times, and 
hence improves sensitivity by a factor which can be estimated by , , \ ^ To 

^ vpulsar duty cycle 

date, VLBI pulsar gating has always used incoherent dedispersion, as the small amount 
of smearing incurred compared to coherent dedispersion has a negligible impact on the 
recovered signal to noise ratio. Pulsar gating is discussed in more detail in Section [4.3.41 
Despite the impressive resolution obtainable with VLBI (< Imas), the small physical 
size of the pulsar emission region means that even the nearest pulsars are completely 



® In this instance, the authors modeled all other contributions and used P{^^'' to measure the distance 
to the pulsar 
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unresolved on VLBI baselines (assuming a pulsar emission height of 1000 km - most likely 
a gross overestimate - yields an angular size of 0.07 /ias at 100 pc). Potential exceptions 
to this rule are the interactions of pulsars which their surrounding environment, such as 
pulsar wind nebulae ( PWN), or interactions with companions, such as PSR B1259-63 



(Johnston et al 



19991 ). although no such detections have yet been published. It is also 



worth noting that VLBI observations of pulsars can be used to obtain very high resolution 
speckle images of scattering disks in the ISM. The first example of such an observation 
(which used the DiFX software correlator) is briefly discussed in Section [4.6.31 

Thus, the main application of VLBI observations of pulsars is to obtain astrometric 
information, either for the purpo se of making a ssociations with other observed structures 
such as supernova remnants (e.g. iBriskenl . |2005| ). or to improve the accuracy of kinematic 
and distance information for use in timing, luminosity, and Galactic electron distribution 
models. It is these latter applications which are the focus of this thesis, and they are 
explored in detail in Chapter [6l 



3 

RADIO INTERFEROMETRY 



3.1 Conceptual overview 

Radio interferometry makes use of the spatial separation of two or more antennas to 
obtain information about smaller angular structures in the radio sky than can be gleaned 
from single-dish studies. Whilst one can probe smaller angular scales with a single dish 
by increasing the dish diameter, this has two undesirable side effects. Firstly, the smaller 
beam means that the survey speed of the instrument is not improved, despite the improved 
sensitivity. Secondly, the cost and technical difficulty imposed by the larger physical 
diameter of the antenna rapidly become prohibitive. The inverse relationship between 
dish size and "primary beam" - the full-width at which the antenna response has dropped 
to one half its peak - is illustrated in Figure ISTTI fa). Figure ISTT b) shows how a pair of 
antennas - the classic two-element interferometer - can discriminate between structure 
which lie within the primary beams of the individual elements. With a full mathematical 
treatment of interferometry deferred until Section [3.51 it is sufficient at this point to note 
the potential degeneracies of a single measurement with the two-element interferometer, 
and to observe that different projected antenna spacings - baselines - would be required 
to determine the structure of the source being observed. 



3.2 Historical development 



Radio astronomy was founded in 1933, when K?'^^ Ja nsky published the detection of the 
Galactic background at low frequency (jjanskvl . Il933l ) . Considerable progress in the new 
wavelength regime, however, was delayed until the end of the second world war, at which 
time large quantities of military radio equipment began to be used for radio astronomy. 
The first interferometric observations were made around this time in Australia, using a 
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Figure 3.1 Resolution of a single dish compared to a two-element interferometer, (a) 
Response of a single antenna element. As diameter D increases, the radiation collected 
at the edges of the dish becomes further out of phase for a given angular offset (j). As 
wavelength A decreases, a fixed time/distance offset at the edge of the dish corresponds 
to a greater amount of phase. The full width half maximum (FWHM) of the antenna 
response - the "primary beam" - is given by \.22\/D. (b) Response of a two-element 
interferometer. The FWHM of the "synthesised beam" is given by the similar expression 
1.22A/i?, where B is the projected distance between the antennas, in a plane perpendicular 
to the observation direction. 
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single receiving element mounted on a sea cliff (jPawsey et alJ . 1 19461 ) . This arrangement 
made use of the path delay provided by the re flection off the sea surfac e; however, arrays 
of separate receiving elements soon appeared (jRvle and Vonberd . 119461 ) . 

These early interferometers measured only the changing intensity of the summed inter- 
ferometer signal - a direct analogue to the optical two-slit experimen t. A major advance 
came with the advent of phase-switching interferometers (jRyld . Il952l ) , which introduced 
a periodic phase inversion to one of the interferometer elements. This allowed the mea- 
surement of the multiplicative term between the elements without the addition of the 
individual squared signals, considerably improving interferometer sensitivity. This, in 
turn, was made redundant through the improved stability of frequency standards which 
allowed direct multiplication of the signals from interferometer elements. 

The ongoing rapid improvements in the capabilities of digital electronic equipment has 
allowed correspondingly rapid improvements in interferometer capabilties. The cost of new 
instruments is now generally dominated by the structural components of the antennas and 
associated infrastructure, making the upgrade of existing instruments with new electronic 
components an attractive proposition. Such upgra des are underway or recently completed 
for the Very Large Array (the Expanded VLA: 
arrays, which are discussed below. 



Perley et al. 



200J) and several VLBI 



3.3 VLBI 



As radio interferometers developed, a natural tendency was to increase baseline length to 
achieve better angular resolution. This trend quickly reached the limits at which informa- 
tion could be could be distributed to and received from antennas in real time with existing 
technology. To overcome this limitation, disconnected stations were equipped with record- 
ing media to store baseband data until it could be brought to a common location, and made 
use of independent freque ncy standar d s. Ea rly ex amples of sci e nce u ndertaken with such 
arrays were published by IClark et al.l (119671 ) and iMoran et al.l (jl967l ) . The technological 



irk et al.l (119671 ) and iMoran et 
limitations on VLBI were progressively lifted, and a number of ad-hoc and part-time VLBI 
arrays functioned around the world from the 1970s onwards. These included the Network 
Users Group in North America^, the European VLBI Network (EVN^, the Asia-Pacific 
Telescope (APTjfl, and the Australian Long Baseline Array (LBAjf|. A purpose-built, 
full-time VLBI array was commissioned in the US in the 1990s - the Very Long Baseline 



^discussed along with other early North American efforts by iKellermann and CohenI (| 19881 ) 
^http:/ /www. evlbi.org/ 
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^http:/ /www.atnf.csiro.au/vlbi/ 



30 



Chapter 3. RADIO INTERFEROMETRY 



Array (VLBAjfl. Throughout the development of VLBI, the use of independent frequency 
standards and non-real time correlation has remained the defining distinction between 
VLBI and regular "connected element" interferometry, although recently pseudo-realtime 
correlation has become possible on some VLBI arrays through the use of high bandwidth 
fibre optic links ( "eVLBI" ) , which is discussed further in Section 14.6.11 For a compre- 
hensive overview of the development of VLBI and radio interferometry in general, see 



Kellermann and Moran (2001 



VLBI allows the highest angular resolution available to imaging in astronomy at any 
waveband, which allows the detailed study of the most distant objects in the Universe, 
as well as objects with small physical extent in the more nearby Universe. Historically, 
targets of interest have included active galactic nuclei (AGN), radio galaxies, masers and 
pulsars. Studying the emission of these objects on the smallest scales has given unique 
insight into the physics that power th em, as well as throw ing up surprising discoveries 
such as apparent superluminal motion (jWhitney et al.l . Il97ll ). One of the most important 



applications of VLBI has little to do with the sources themselves, however, and is the cre- 
ation and maintenance of a stable, quasi-inertial reference frame upon which astronomical 
positions can be based. Historically the domain of opti cal astronomy (and defined most 
recently by the FK5 reference frame: iFricke et al.l . Il988l ). this responsibility was assumed 



by VLBI when the International Celestial Reference Frame (ICRF) superse ded the FK5 



reference frame in 1998. Defined by 212 distant radio sources (iMa et al 



19981 ). it has been 



extended to include several hundred additional "candidate" and "new" sources |Fev et al 



20041 ). VLBI measurements also form an important contribution to the International Ter- 



restrial Reference Frame (ITRF), which along along with the ICRF is discussed further in 
Section 13.81 



3.4 Instrumentation and hardware 

The purpose of an element in a modern radio interferometer is to obtain a digitised time- 
domain representation of the radio energy within a desired frequency band incident upon 
the element, adding as little noise as possible in the process. Once this has been achieved, 
digital signal processing can be used to compute the correlation between array elements, 
which along with a knowledge of the array geometry can be used to generate a model of 
the observed sky in the first the spatial frequency domain, then the image domain. This 
section describes the hardware required - a mathematical description of the foundations 
of interferometry is presented is Section 13.51 The description below assumes a typical 
^http://www. vlba.nrao.edu/ 
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Figure 3.2 Hardware components of an interferometer. 



reflecting concentrator (dish) type radio telescope, referred to throughout as an antenna. 
Figure [3^2] shows the components discussed below. 

The electric field present at the focal point of an antenna is converted to an electrical 
voltage by a sensor referred to as a feed. Typical feeds consist of two components sensitive 
to two (usually linear, although circular is possible) orthogonal polarisations. The feed 
itself typically resides within a feed horn, whose purpose is to shape the illumination of 
the feed to the antenna surface and hence optimise the antenna primary beam response. 
The time-varying voltage signal v{t) at the feed is linearly proportional to the electric 
field E{t) present at the feed. 

The voltage signal induced in the feed is extremely small, as the radio power influencing 
the voltage is itself extremely small. Thus, the first stage of signal manipulation is to 
amplify the signal using a Low Noise Amplifier (LNA). This first stage of amplification 
dominates the noise budget of the receiver system, and so considerable effort is made to 
minimise the noise contribution of the LNA. The dominant noise contribution comes from 
thermal electron motion, and so feeds and LNAs are often housed in a cooled chamber - 
a dewar - as shown in Figure 13.21 
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The calibration of received power is generally undertaken through the injection of 
a known noise source prior at the beginning of the receiver system, prior to the LNA. 
By modulating the noise source on and off, the rise in measured power can be used to 
calculate the total power. This procedure is particularly important for VLBI, where there 
are no unresolved, constant flux sources that can be used for absolute flux calibration, as 
discussed in Section 13.61 The noise source is usually a thermally controlled resistive load, 
which is coupled to the input of the LNA. 

Once the received signal has been amplified, it must generally be downconverted from 
the sky frequency to a lower frequency where digital electronics can be efficiently used. 
This is achieved by a series of mixing operations, often with filtering applied at each stage 
to prevent aliasing of the signal. The reference oscillator for the mixing stages is generally 
based on an atomic oscillator, such as rubidium or hydrogen (commonly used for VLBI due 
to the excellent frequency stability). This "local oscillator" is distributed to all antennas 
in a connected element array, but this is not feasible for VLBI arrays. 

Once the data is at or close to "baseband" (frequency range starting at Hz), it can 
be sampled and digitised, at which point digital signal processing can be applied. The 
sampler /digitiser (and recorder, for VLBI) is collectively referred to as a "digital backend" , 



and modern examples incorporate features such as digital filtering (e.g. Ilguchi et al.l . l2005l ) 
to improve the data quality. The digitised data is then transported to the correlator, 
which produces the sampled visibilities (scaled fractional correlation between antennas) 
as a function of frequency and time. The correlator is the focal point of the hardware 
chain in terms of this thesis, and correlators in general and the DiFX software correlator 
in particular are discussed in more detail in Chapter [H 

3.5 Mathematical description 

All of the salient concepts of interferometric Fourier synthesis imaging can be succinctly 
illustrated using a two-element interferometer, so this simplified example will be presented 
here. It should be noted t hat this sec t ion d raws heavily on the explanation of interfer- 
ometer theory presented in bhompsoJ ligg^). Throughout, the standard interferometric 



assumptions of a far-field, spatially incoherent source are made. A diagram of a two- 
element interferometer is shown in Figure [3l3l The elements are separated by a baseline 
b, and are both pointing in the direction of the unit vector s. The geometric time delay 
between the signal arriving at antenna A and antenna B is given by Tg{t) = b • s/c, where 
c is the speed of light. It is shown as a function of time since, depending on the choice of 
reference axis, either b or s will change with time as the Earth rotates. 



3.5. Mathematical description 
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Considering the astronomical signal initially to be a monochromatic wave of frequency 
f considerably simplifies the analysis, reducing the correlator function to a simple multiply 
and accumulat^, while still demonstrating the Fourier relation between interferometer 
observables and the actual sky brightness distribution. In this case, the response of a 
single element will be of the form V = v cos 2TTvt, and so the correlator output can (in the 
absence of any delay tracking) be written as: 



r{t) = < VA{t)VB{t) > (3.1) 
= vavb cos2T:v{t — Tg) cos2t:v 
= vavb cos ATTvi — 2-1: vTg + vavb cos 27ri/rg 
— VAVB cos2TruTg 

Thus, assuming the averaging time is long compared to the term at twice sky frequency, 
but short compared to the changing Tg, the correlator output will be a sinusoid-like func- 
tion (since Tg is not varying linearly, it is not a true sinusoid) with amplitude VaVb, which 
is proportional to the received power. 

If radio brightness is represented by /, then the radio brightness in the direction of s 
can be written as /(§). Brightness, which is the desired quantity when mapping the radio 



the non monochromatic case will be considered in the following chapter 
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sky, is measured in power per unit area, per unit bandwidth, per unit solid angle, and can 
be placed in units of W Hz~^ sr~^. From this it follows that from a source element 
dVt in the direction s, an antenna with effective area ^(s) accepting bandwidth Az^ will 
receive power equal to A{s)I{s)^v d^. Since the correlator output was proportional to 
received power, by neglecting constant gain factors the correlator output from solid angle 
dVl can be written as: 

dr = ^(s)/(s) Az^ dVL cos l-KUTg (3.2) 
The correlator output can then be written as integration over the whole sky (5 = 47r 

r = Ai^ [ A(s)/(s)cos '^^^^ " ^ dn (3.3) 
Js c 



steradiansH: 



In general, interferometric images are made in a relatively small solid angle, as con- 
strained by the antenna primary beams. Thus, it is convenient to rewrite s as s = Sq + cr, 
where the image centre sq is fixed and is referred to as the phase centre. Substituting this 
into equation 13.31 yields: 

r = Ai/cos / ^ cr /(o-) cos dO. 3.4) 

c Js c 



. iTTvh - So f ^r/ \ ■ 27ri/b • cr 

Az^sm / A{a)l{a)s\n dil 

c Js c 



At this point, it is necessary to introduce the term visibility. The visibility is a measure 
of the coherence of the radio sky brightness distribution, modified by the antenna response, 
between an antenna pair. It is a complex quantity defined as: 

V = \V\e''f'^ = 4- / A(cr)/(o-)e2™*'-'^/^dJ7 (3.5) 
^0 Js 

where Aq is the antenna response at the phase centre. Subsequently, it will be shown 
that the visibility is, under certain assumptions, actually the Fourier transform of the radio 
sky brightness /(cr). First, however, the relationship of correlator output to visibility will 
be shown. By separating the real and imaginary components of Equation l3.5^ the following 
expressions are obtained: 

Aq\V\cos(I)v = [ A{(t) I (a) cos ^'''''^''^ dn (3.6) 
Js c 



'In practice, however, the hmitations of the primary beam of the antenna elements means that A{s) is 
only non-zero for a small solid angle 
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and 

Aq\V\s\.u4)v = - [ A{(t) I {a) sin ^^^^^ " dn (3.7) 
Js c 

which can then be substituted into Equation 13.51 to obtain: 

r = yloAHF|cos(^^^^^^^-<^y) (3.8) 

Thus, a correlator with no delay tracking measures the visibility, modulated by a 
fringe pattern with frequency dependent on frequency and baseline, and an amplitude 
dependent on the antenna response. The removal of this fringe pattern by a realistic 
correlator to enable direct measurement of visibility, as well as dealing with the effects of 
finite bandwidth and frequency down-conversion, are discussed in Chapter HI 

Once the complex visibility has been sampled by an interferometer, the desired result 
is generally an image of the radio sky brightness distribution. In order to show the validity 
of synthesis imaging, it is necessary to define the coordinate systems for both visibility 
and sky brightness. These are illustrated in Figure [331 The three orthogonal coordinate 
axes for a visibility point are (u,v,w), with the corresponding axes in the brightness 
distribution being (l,m, n). The natural units for the visibility axes are wavelengths cor- 
responding to the observed frequency. This convention will be followed for the remainder 
of this thesis. In general radio astronomy convention, the w axis points in the direction 
So, while u and v are orthogonal to each other and the w axis, and are on the plane 
containing the w axis and celestial East and North respectively. The synthesised image 
is assumed to be two-dimensional, and so positions on the sky are described in terms of I 
and 771. The vector quantity s can be expressed as II + mm + nh. 

The quantities b, s, sq and Q referred to in the above derivation can then be expressed 

as: 

i^b • s 

= ul + vm + wn (3.9) 

= w (3.10) 

c 

,^ dl dm dl dm 

dn = = , 3.11 

Substituting Equations I3.9l[3m into Equation 13.51 yields: 

V{u,v,w) = ^ r r A(;,m)/(Z,m)e-2'^^["'+'^"^+"'(^^^^^-i)l^iS^ 

(3.12) 

Since, as components of a unit vector, + m^ + = 1, the integrand is necessarily 
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m 




*- East 



Figure 3.4 Definition of interferometer coordinate systems used throughout this thesis. 
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where + m? > 1. 

To obtain the radio sky brightness distribution 1(1, m), it is necessary to invert Equa- 
tion I3.12[ Whilst possible in the generalised form given above, it is conceptually difficult 
and computationally expensive. If the visibility equation given above can be reduced to 
the form of a two-dimensional Fourier transform (necessitating t he removal of the w term ) , 
the inversion can be carried out using a Fast Fourier Transform (jCoolev and Tukevl . ll965l ). 
which is much more computationally tractable. If / and m are constrained to be small - 
restricting the field of view to a small solid angle - then the term \/l — P — rm? reduces 
to 1, eliminating w from Equation 13.121 and leaving: 



V{u,v) 



1 



A{1, m)I{l, m)e-2^*["'+^™J dldm 



(3.13) 



Equation 13.131 can be inverted via the Fourier relation to yield: 
1 



-A{l,m)I{l,m) 



(3.14) 



Thus, the correlator output can be transformed, under certain assumptions, to yield 



the original sky brightness. For a more detailed descripti on of the 



small-field assumption and methods of circumventing it, see iThompsonI (jl999l ). 



imitations of the 



3.6 Calibration and editing 



Before the transformation of visibilities into a sky brightness distribution can take place, 
incorrect and uncorrectable visibilities must be removed (a process known as flagging), 
and the remaining visibilities correctly calibrated. Visibilities which must be flagged are 
typically those which have been affected by equipment failure, radio frequency interference 
(RFI) or severe propagation problems (rapid ionospheric visibility at low frequencies, water 
vapour variations at high frequencies, and physical problems such as antenna shadowing). 
Flagging can be automated (based on logs of telescope and electronic failures, for example, 
or on the visibility values themselves) or manual, where visibilities are inspected by eye 
and questionable samples removed . Flagging is an inexact proced ure which has been the 



200e 



Ekeri 



19991). 



subject of numerous analyses (e.g. iMiddelberg . 

Calibrating the interferometer output amplitude to obtain the true visibility amplitude 
requires, as shown in Equation 13. 81 a knowledge of the peak antenna response. As discussed 
in Section 13.41 this can be achieved through the use of noise calibration at the individual 
antennas. For short-baseline arrays, this calibration is typically reflned by including an 
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observation of a known "flux calibrator", a compact, unconfused source known to possess 
stable flux. Antenna-based correction factors can be derived by comparing the measured 
visibilities with the known values. This is a luxury not available to VLBI arra ys, since 



there are no constant flux sources which are still compact on VLBI baselines (IWalkerl . 



19991 ). Thus, accurate logging of antenna noise calibration is especially crucial for the 
absolute, as well as relative, calibration of VLBI observations. 

At this stage, it is also possible to make corrections to the visibility phases, to compen- 
sate for dispersive and non-dispersive delays which were not correctly accounted for in the 
correlator model - that is, the first term in the cosine of Equation 13.81 was not correctly 
subtracted. Typically, this involves measured atmospheric and ionospheric quantities that 
were not available at the time of correlation, and are most commonly required in VLBI ob- 
servations. Other potentially correctable errors which can occur during VLBI correlation 
include geometric model errors due to incorrect Earth Orientation Parameters (EOPs). 
These values, which collectively describe the orientation and rotational phase of the Earth 
to very high precision, are made available by the International Earth Rotation Service 
(lERqfl), and are calculated using intensive geodetic VLBI observations (see Section [3.81 
below) . 

Once calculated, the application of corrections is simple, amounting to a simple scaling 
of visibility amplitudes or rotation of phase. Typically, the rate of change of corrections is 
small, and so the error caused by the finite visibility sampling duration is not a concern. 
For very rapidly varying corrections, or for errors which cannot be modeled and subtracted, 
however, the sampled visibilities remain in error. Residual visibility errors add noise to 
an inter ferometric image, but may be corrected (in some circumstances) in the imaging 
stage as discussed below. 

3.7 Imaging 

This section is meant as an overview of the main techniques used in VLBI to produce radio 
images. It does not cover complicating factors such as the (potentially inhomogeneous) 
antenna primary beam responses, wide-field effects, the effect of spectral and temporal 
averaging, or gridding techniques. Fo r a c omplete reviews of radio synth esis imaging, see 



(for example) [Thompson et al.l (|1994| ) and iPearson and ReadheadI (|1984| ) 



Once an edited, calibrated visibility set has been produced, it may be transformed into 
an image using the Fourier relation shown in equation 13. 141 Generally speaking, however, 
the sampling of visibilities in the uv plane will be incomplete (especially for VLBI, where 



' ht tp : / /www . iers . org/ 
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small numbers of antennnas often mean extremely sparse sampling), which introduces 
considerable distortion in the transformed image. Essentially, the measured visibilities 
can be considered to be the product of the true Fourier transform of sky brightness with 
a sampling function S{u,v), which is 1 (or indeed any non-zero weight, if weighting is 
applied in the uv plane) where visibilities were measured and elsewhere. By virtue of 
the Fourier transform duality between multiplication and convolution, this leads to the 
transformed image effectively becoming the convolution of the true sky brightness with the 
Fourier transform of the sampling function ^ {S{u, v)} = s{l, m). This directly generated 
map is generally referred to as the "dirty" map, and the Fourier transform of the sampling 
function, s{l,m) is known as the synthesized beam. 

Thus, to obtain the true sky brightness distribution, it is necessary to "deconvolve" 
the effects of limited visibility sampling. In general, this takes the following iterative form: 

1. Adjust sky brightness model (add/subtract /move components). 

2. Compute visibilities expected due to this model. 

3. Subtract expected visibilities from observed visibilities to obtain a "residual" bright- 
ness distribution. 



Once the best source model has been obtained, it is convolved with a "restoring" 
beam, which is typically taken as an elliptical gaussian truncated at the first null of the 
synthesized beam, and the remaining residual map added. This generates a "clean" map 
with similar resolution to the original dirty map, but with the sidelobes generated by the 
synthesised beam minimised. 

Various procedures have been developed to make use of a priori information to help 
guide step 1 above, and converge towards the best (and simplest) global solution for the 
true sky brightness, which is complicated by the fact the the measured visib ilities are 



corru pted by noise. Popular exa niples of imagin g algorithms include CLE AN (IHogboinl . 



19741 ) ■ least squares model fitting (jPearsonl . Il999l ). maximum entropy (MEM; 



Abies 



19741 ^ 



and various derivatives and combinations of the above. In this thesis, both model fitting 
and CLEAN were used. 

During the imaging process, it is possible to make use of additional relationships be- 
tween groups of visibilities to improve the a priori calibration and remove errors in the 
visibilities. These "closure" quantities result from the relationship between the true and 
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observed visibilities, which can be written as shown in Section 13.51 



K^n = 9m9lVmn (3.15) 

ICnle^"^™ = |ffm||ffn||K.n|e'(^«--'^-+^-) (3.16) 

where gm = |grrt| exp i (/i)^^ is the complex gain associated with antenna m. Considering 
only terms in the exponent of Equation 13. 161 vields: 



I obs. 



mn 



+ (pmn (3.17) 



and considering the summation of the phase terms from the common baselines of three 
antennas m, n and p gives the phase closure relationship for the "closure phase" 4>rnnp- 

0^n| = + + = ^mn + 4>nv + 4'vm (3.18) 

Equation 13.181 shows that the closure phase is a quantity which is independent of the 
antenna-based gains, which include the effect of varying atmospheric propagation and 
other antenna-based errors, and as solely determined by the true visibility phase. Thus, 
if a sufficiently accurate model of the brightness distribution is available, deviations in 
phase from the model phase can be attributed to antenna-based errors and adjustments 
made to the value of for each antenna to minimise the difference between model and 
observed visibility phases. This procedure is known as "phase self-calibration". 

By considering the magnitude terms of Equation l3.16l the following equation for closure 
amplitude iV^mnpgl can be constructed from baselines of four antennas m, n, p and q: 



iTAclosei _ I ''ran 1 1 pq I _ | '"mnM ^ pq\ /„ -,„s 

Equation 13.191 shows that the closure amplitude, like closure phase, is independent of 
the antenna-based complex gains, and similarly can be used to compute a correction to the 
gains to minimise the difference between model and observed visibilities. Unsurprisingly, 
this procedure is known as amplitude self-calibration. 

While self-calibration is a useful tool for removing antenna-based errors and improving 
the fidelity of radio synthesis images (usually measured by dynamic range - the ratio 
of peak image flux to image noise), its principal drawback is that it can only attempt 
to reduce the difference between the observed visibilities and the model visibilities. If 
the input model does not resemble the actual sky brightness distribution, application of 
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self calibration will not improve image quality, and may limit the potential to improve 
the model in the future. Furthermore, if the closure quantities are extremely noisy, self 
calibration may act to fit noise to the model. For more details on the use of self-calibration 
in interferometric imaging, see lCornwell and FomalontI (|l999l ). 



3.8 Astrometry and geodesy 



As shown in Equation the uncorrected phase term of correlator output contains the 
true visibility phase for a source, corrupted by a geometric term. Actual correlator im- 
plementations, as discussed in the following chapter, attempt to remove this non-intrinsic 
term as precisely as possible. However, the ability to do this depends on the exact knowl- 
edge of both the b and s vectors, and thus on the exact location of the source and antennas. 
Errors in the assumed values for b or s will cause the correlator output phase to differ 
from the true visibility phase. 

The use of the visibility phase error to calculate error in assumed source position (hold- 
ing the baseline vector fixed) is known as astrometry, while calculating errors in assumed 
baseline vectors is known as geodesy. An excellent ov erview of t he th eory and practical 
difficulties of VLBI astrometry and geodesy is given in iFomalontI (jl999l ) . It should be im- 
mediately apparent that disentangling the phase contributions due to unknown baseline 
and source vector errors simultaneously is a challenging task, and multiple baselines and 
sources are required to make such a global solution. Further complicating the problem 
are the residual visibility phase errors due to the unmodeled atmosphere and ionosphere. 
The contribution from the ionosphere is frequency-dependent and can be estimated and 
subtracted if observing frequencies span a wide fractional bandwidth and sufficiently high 
signal to noise ratios can be obtained, and astrometric/geodetic observations often use 
a dual-band "S/X" (2.3 GHz and 8.4 GHz ) receiver for the purpose of obtaining widely 
separated bands (see e.g. 



Petrov et al. 



20081 ) . The contribution from the atmosphere (pre- 



dominantly the wet and dry troposphere) is independent of frequency, however, and so 
improvements to the a priori model can only be r nade through the use of of station weather 
information, such as water vapour radiometers (|Roy et al.l . 120061 ). 

The principal difficulty confronting astrometric and geodetic observations is that the 
powerful self-calibration technique cannot be used freely - since it acts to converge visi- 
bilities with a pre-existing model, incorrect usage will corrupt the measurements of source 
and antenna positions which are being sought. Sources used for astrometry and geodesy 
should ideally have no detectable structure, simplifying the prob lem, but suitably bright 



examples of such sources are rare (see e.g. 



Gontier et al 



2001 



). Thus, correct calibra- 



42 



Chapter 3. RADIO INTERFEROMETRY 



tion of visibility phases through the careful modeling of propagation effects and source 
structure is of the utmost importance in astrometric and geodetic observations. 

The fundamental reference frames which are the product of absolute astrometry and 
geodesy are the ICRF and ITRF, introduced in Section 13.31 They are maintained and 
improved by the lERS. The ongoing maintenance of these frames allows other users to 
undertake relative astrometry and (less commonly) geodesy, under the assumption that 
positions of the defining members of these frames are correct. This procedure of relative 
VLBI astrometry forms the basis of the measurement of pulsar position measurements 
made in this thesis, and as such is discussed in further detail in Chapter [5l Geodesy is 
not a fundamental component of this thesis, but will be mentioned where appropriate as 
it impacts on the astrometric measurements. 

When undertaking relative, as opposed to absolute astrometry, one of the major prac- 
tical differences is the use of single-band phase, as opposed to multi-band delay which 
is commonly used for absolute astrometry and geodesy. Even the best available a priori 
delay modeling is uncertain at the level of nanoseconds, which corresponds to multiple 
turns of phase at typical observing frequencies. Thus, for absolute astrometry the de- 
lay is generally estimated by the phase gradient between observing bands - the so-called 



Fomalont 



"mult i-band delay". The concept of multi-band delay is explained in detail in 
(jl999l ). When conducting relative astrometry, however, the delays are calibrated to higher 
precision through observations of a known source, meaning the single-band phases no 
longer suffer from wrap ambiguities. This allows a higher precision estimation of relative 
position. However, the accuracy of the delay solution transfer from calibrator to target in 
relative astrometry depends on the angular separation of the two sources, and the gradient 
in residual delay between them. Clearly, it is desirable to minimise this angular separation, 
in addition to attempting to obtai n the best possible modeling to minimise the residual 



error gradient. iPradel et al.l (j2006l ) present an comprehensive simulation-based analysis 
of errors in relative astrometry, including the effect of different error sources such as the 
wet and dry troposphere at different angular separations. 



4 

DIFX: AN FX STYLE SOFTWARE 

CORRELATOR 



An in -depth description and analysis of the DiFX code was pubhshed in iDeller et al 
(|2007l ). and the remainder of this chapter draws upon the n iaterial p r esent ed there, as 
well as the description of correlator functionality presented in iRomnevI (|l999l ). 



4.1 Mathematical description of a correlator 

Chapter [3] showed that a radio interferometer measures the spatial frequency components 
of the observed sky brightness distribution, and that this information, under certain re- 
strictions, could be transformed into an estimation of the actual sky brightness distribu- 
tion. However, this treatment was valid only for the unphysical case of a monochromatic 
radio source, and neglected the effect of data sampling, which is crucial to the operation 
of modern interferometers. In this section, I show the effects of relaxing these restrictions, 
and the functionality which must be incorporated into a realistic correlator to compensate. 
The description will be appropriate for an unconnected interferometer, and simplifications 
that can be adopted for connected interferometers will be noted where applicable. The 
mathematical derivation will be appropriate for an FX-type correlator (see Section [4.2.11 
for the definition of an FX-type correlator), and the advantages and limitations compared 
to an XF-type correlator will be highlighted. 

4.1.1 Quasi— monochromatic formalism 

It is convenient to represent the electromagnetic radiation from the source as a quasi- 
monochromatic plane wave incident on the antennas, described (for a given antenna p) 
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by: 

Ep{t + Tp) = A{t + rp)e^2.'^o(*+^^) (4.1) 

The electric field component is given by the real component of Equation 14. 11 It is clear 
that for a realistic, band-limited process A{t + Tp), the electric field component will be 
a band-limited signal, centred on i^q. A{t + Tp) is a complex signal which must also be 
covariance-stationary, ergodic and stochastic. The term Tp s. 
represents the delay between antenna p and a reference point ^. 



lown in Equation 14.11 above 



Since A{t + Tp) has been asserted to be band-limited, it can be written as: 

/oo 
s(^i,)e^^Mt+rr,)^^ (4.2) 



00 



with the proviso that s(j^) = outside some limiting band |z^| < B. 

As shown in Section [3.51 tlis correlator's function is to compute the expectation of the 
cross correlation of the two received signals from antennas p and q: 

rp,{T) = {Ep{t + Tp)E;{t + T,)) (4.3) 

= {Ep{t + Tg)E;{t))) 

= {A{t + Tg)A*{t))e'^'""'^' 

= lp,{r)e'^^''°^^ 



Throughout this thesis, the superscript * will be taken to mean complex conjugatio: 
Since A{t) is ergodic, the expectation can be approximated by a time average. As in 
Chapter [3l Tg represents the geometric delay Tp — Tq between antennas p and q. The 
term ^pq{T) is the covariance function of the stochastic process A[t), and represents the 
unmodulated correlation between the electric fields on this baseline - the visibility, as 
shown in the previous chapter. This visibility is modulated by the time-varying phasor 

Ultimately, the visibility is required as a function of frequency, which can be obtained 



^Usually the geocentre for VLBI, or some point near the centre of the array for connected element 
interferometers 

^The real (non-complex) correlator described in Chapter |3] did not require the second signal to be 
conjugated in order to calculate the correlation function, for the obvious reason that all signals were real 
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by substituting the spectral representation for A{t) into the expression for ^pgir) as shown: 

7p,(t) = {A{t + Tg)A*{t)) (4.4) 
= { / Sp{u)sl{u')e'^''^^''-'''^'+''^oUi^diy') 



— oo J — oo 
oo 



Spg{i^)e'^''''^^diy 



oo 



where Spq{v) = Sp{i^)s*{i') is the cross-power spectrum of A{t), and can now be shown to 
be simply the Fourier transform of the covariance function: 



Spgii^) = / 7pg{T)e-^''"'^dT (4.5) 

J — oo 

Thus, in order to obtain the term Spq{i'), which is the desired visibility term, a realistic 
correlator must convert the time domain signal to the frequency domain and account for 
the changing geometric delay Tg and other practicalities such as frequency downconver- 
sion and sampling - in effect, compensating for the phasor term e^^Tri^oTg^ These steps are 
detailed in the following sections, assuming station-based corrections. Station-based re- 
moval of the phasor fringe term is generally only computationally feasible in an FX-style 
correlator architecture, as explained in Section [4.2.11 



4.1.2 Baseband conversion and sampling 

Despite continual improvements in digital electronics which have produced cheaply avail- 
able, high-speed (multi-Gsample s~^) digital samplers, some form of frequency down- 
conversion is ubiquitous in all arrays, with the exception of dedicated low-frequency in- 
struments. For observations at wavelengths of several cm and shorter, this situation is 
unlikely to change in the near future. Typically, frequency conversion is a multi-stage 
process, utilising broad front-end filters, sharp intermediate-frequency (IF) filters and 
potentially digital filtering after sampling. The mathematical treatment below assumes a 
single conversion step directly to baseband, which is unrealistic but ma thematically con 



venien t. For an in-depth discussion of frequency conversion systems, see [Thompson et al 



(I1994j) 



The application of a real, single sideband mixer with frequency i^q to the received signal 
E(t) from Equation 14.11 yields the following signal: 

V'{t + Tp) = Ep{t + Tp) X cos 2TiUQ = ^ Li2nuoT, ^ ^i2nM^t+T,)] (^^g) 
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The application of a low-pass filter removes the high-frequency term from Equa- 
tion 14.61 At this point, it should be noted that in a connected element interferometer 
with a common clock, the phase of the downconversion signal can be varied continuously 
to compensate for the changing geometric phase (effectively replacing the oscillator signal 
gi27rj/oi -^itJi ^i2nuoit~Tp)y This renders the fringe-rotation and fractional sample correc- 
tions detailed in Section 14.1.31 unnecessary. However, as this approach is not presently 
implemented on major VLBI arrays, the following analysis will assume that no phase 
compensation has been performed at the downconversion stage. 

After frequency downconversion, the resultant bandlimited, baseband signal Vp{t) is 
sampled. Since the process A{t) is bandlimited with bandwidth B, the sampling interval 
must be no longer than At = 1/2B in order to prevent aliasing of the sampled spectrum. 
Quantisation can be as coarse as 1 bit, although 2 bits has been widely used and newer 
instruments have tended to push towards higher precision in order to mitigate the effects 
of RFI, and obtain higher spectral dynamic ranges for spectral line observations. At the 
low bit precision typically used, the sampling function S is typically non-linear, allowing 
an optimum placement of levels given the Gaussian distribution of V. The sampled signal 
Rp{n) can then be written as: 



For the sake of clarity, the following discussion of geometric compensation will use the 
continuous notation for the baseband time series Vt, and the limitations placed by the 
sampling process will be highlighted where appropriate. 

4.1.3 Geometric compensation and channelisation 

In order to form the correct correlation between the array antennas, the geometric delay 
Tp for each antenna must be removed. Initially, this is implemented as an integer-sample 
shift of the sampled signal Rp{n) by N samples, where N is the rounded integer value of 
Tp/At. The shifted baseband signal is thus Vp(t + Tp - nAt) = A{t + Tp - nAt)e^'^^'"^''p . 
Since the integer-sample shift cannot exactly compensate for the geometric delay r^, a 
"fractional-sample" error Cp = Tp — nAt remains at this point. A diagram of the sampled 
signal Rp{n) which illustrates the geometric delay Tp, the integer-sample delay NAt and 
the fractional-sample error for an antenna p at a given time is shown in Figure 14.11 

The compensation of the fringe phasor e*^'^''"'^'' can now performed. This step is gen- 
erally referred to as fringe rotation. Conceptually this is simply applied as a complex 




(4.7) 
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Figure 4.1 Geometric delay for a sampled signal Rp{n) for an antenna p at a given instant 
of time. The exact geometric delay Tp is shown in blue. The nearest integer-sample 
delay (shown in green) can be compensated by shifting the sampled datastream, leaving 
a fractional sample error Cp (shown in red). 



multiplication to Vp by the function ^-^'^'^'^o^p ^ yielding the fringe-rotated, shifted base- 
band signal Vp(t + ep)e~'^'^'^'^°'^p = A{t + gp), which is complex- valued. An alternative to a 
complex multiplication is to implement the multiplication with separate real-valued cosine 
and sine components. This ensures that the signal remains real-valued and allows the use 
of real-valued Fourier transformations. The "cosine" and "sine" arms of the correlator 
are then combined after Fourier transformation by applying a Hilbert transform to the 
sine arm befo re addition. For a n exa mple of this form of implementation of a complex 
correlator, see [Thompson et al.l (j 19941 ). 



The shifted, fringe-rotated signal may now be transformed into the frequency domain 
as shown: 

9" {A{t + ep} = s^e'^^"'" (4.8) 



In practice, this step is implemented using a Fast Frequency Transform (FFT: ICoolev and Tukevl . 



19651 ). This imposes a segmentation upon the data, since the FFT windows the time do- 



main data with a window size equal to twice the desired number of spectral points. 

The final geometric compensation step is the removal of the phase gradient generated 
by the presence of fractional sample error. Since the maximum magnitude of e is half 
the sampling time (1/(45) seconds), the phase induced at the upper edge of the band 

= B) can range between — vr to vr. Fractional sample correction is implemented as a 
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Table 4.1. Maximum decorrelation incurred due to "Post-F" fringe rotation 



Observation 


Maximum 


Frequency 


Spectral points per 


Maximum 




baseline (km) 


(MHz) 


16 MHz band 


decorrelation (%) 


LBA low frequency eontinuum 


1400 


1600 


128 


0.003 


LBA high frequency continuum 


1700 


8400 


128 


0.13 


VLBA low frequency continuum 


8600 


1600 


128 


0.12 


VLBA high frequency continuum 


8600 


22200 


128 


21.1 


LBA water masers 


1700 


22200 


1024 


47.6 



complex multiplication of the frequency-domain data. After fractional sample correction, 
the spectral representation of the original signal s{u) has been recovered exactly. 

Finally, it is important to note at this point that realistic correlators approximate 
the true geometric delay function Tp over short time intervals by a polynomial expansion. 
Over a the time duration of a single FFT window (typically microseconds, for bandwidths 
B ~MHz and spectral points n ~hundreds), this can be (and usually is) approximated 
by a linear expansion. If the difference in fringe phase across an FFT window 6(f) = 
2ttuq (Tp(ti) — Tp{ti + 2nAt)) is small, the fringe term may be approximated by a constant 
over the entire FFT window. Since the Fourier transform is a linear transformation, this 
single complex factor then may equivalently be applied after the FFT, rather than to every 
individual sample before the FFT. In practice, this allows the fringe rotation operation 
to be merged with the fractional sample correction. "Post-F" fringe rotation allows a 
considerable saving in computational load, since for every FFT window it removes the 
need for 2n multiplications and 2n trigonometric operations, at the cost of n complex 
additions and one trigonometric operation. However, if the FFT window is long (high 
spectral resolution) and/or the fringe rate is very high (high frequencies and/or long 
baselines), the change in fringe phase over an FFT window may incur an unacceptable 
amount of decorrelation. The decorrelation can be calculated as smc{5(j)/2). Table ST] 
shows the maximum decorrelation which would be incurred through the use of post-F 
fringe rotation for some typical observations. 

4.1.4 Cross— multiplication and accumulation 

After each station has been individually delayed, fringe-rotated and corrected for frac- 
tional sample error, the desired visibility term Spq{v) can be obtained through a complex 
multiplication of the corrected baseband station data, after conjugating the second data 
stream: 

Spq{u) = sp{u)sl{u) (4.9) 
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This is repeated for all baselines, and the result is accumulated for a desired time in- 
terval, typically seconds for narrow-field VLBI. If p = g, the result is an autocorrelation, 
which is necessarily real-valued. At the end of each integration period, the visibility data 
is normalised and stored. The normalisation can be performed in an a priori manner (in 
theory it is dependent only upon the integration time, since the normalised autocorrela- 
tions should tend to unity), but quantisation threshold errors make corrections based on 
the measured autocorrelation necessary. The corrections can be calculated and applied at 
correlation time, or in later processing. 



4.2 Correlator implementations 



4.2.1 FX vs XF correlators 

As shown above, it is possible (and mathematically simplest) to form the visibility out- 
put of a correlator by Fourier transforming the station data, and then cross-multiplying 
the corrected data streams. This approach of channelisation (F) followed by cross- 
multiplication (X) is known as an FX-style correlator. Through the Fourier duality 
between convolution and multiplication, it is easy to see that another possible imple- 
mentation of a correlator would be to form the convolution of the time series for each 
baseline and Fourier transform the result. This is the so-called XF style of correlator. 
Until relatively recently, all correla tors have been i mplem ented as XF correlators - the 



FX concept was first implemented by IChikada et al.l (j 198 71 ) . As the correlator developed 



for this thesis is of the FX type, details of the implement ation of an XF- style correlator 



will not be given here - the interested reader is directed to iRomneyl (119991 ) and references 



therein. However, some of the functional differences between XF and FX correlators are: 

1. In an XF correlator, it is possible (indeed, necessary, to make the design computa- 
tionally viable) to accumulate the lag values and perform a single Fourier transform 
on the accumulated lags before storing the visibilities. 

2. XF correlators can operate at much lower precision than FX correlators and can 
generally make use of optimised, low-precision integer arithmetic. The presence of 
a Fourier transformation before accumulation (and the use of station-based fringe 
rotation) means that an FX-style correlator requires higher precision operations, al- 
though the FX architecture requires less operations overall. The use of low-precision 
operations (including fringe rotation) is generally necessary to make the architecture 
viable. 
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3. Due to the accumulation performed before transformation to the frequency domain, 
the fractional sample error cannot generally be approximated by a constant in an XF 
correlator, which means that correction of fractional sample error in the frequency 
domain in a station-based manner is not possibljf]. The usual alternative employed 
by XF correlators is baseline-based alignment of data, fringe rotation and fractional 
sample correction. This allows alignment of the data streams to within half a sample, 
as opposed to one sample with station based alignment, which when coupled with 
coordinate d fringe rotat ion allows minimal decorrelation due to fractional sample 
errors (see lRomneylll999l for details). Low precision fringe rotation (which cannot be 
employed in station-based fringe rotation due to the possible of spurious correlation 
between harmonics of the fringe function) can also be used in this approach for 
higher efficiency. 



4.2.2 Hardware platforms 

Both XF- and FX-style correlators have traditionally been highly application-specific de- 
vices, based on purpose-built integrated circuits (Application Specific Integrated Circuits 
- ASICs). In the last 20 years, Field Programmable Gate Arrays (FPGAs) have become 
popular components in correlator design s, with one promin ent example being the Very 
Long Baseline Array (VLBA) correlator ([Napier et al.l . 11994 ). FPGAs are reconfigurable 
or reprogrammable devices that offer more flexibility than application-specific integrated 
circuits (ASICs) while still being highly efficient. However, FPGAs still require a greater 
familiarity with the underlying hardware than coding on a general purpose CPU, and 
configurations can be specific to the particular FPGA, limiting portability. Through- 
out the remainder of this thesis, correlators based on ASICs and/or FPGAs will be re- 
ferred to simply as "hardware correlators". Mo dern VLBI hardwar e correlators include 
the NRAO Very Long Baseline Array co rrelator iNapier et al. . 1994 ): the Joint Institute 
for VLBI in Europe (JIVE) correlator (jCassd . Il999l ): the Canadian NRC S2 correlator 



( Carlson et al 



I999I'): the Ja panese VLBI Space Observatory Programme (VSOP) cor- 



relator (Horiuchi et al 



correlator ( Wilson et al 



200d): and the Austraha Telescope National Facility (ATNF) S2 



19961 ^ 



Generally speaking, hardware correlators allow the most efficient design possible, by 
tailoring the hardware to the processing required. However, this efficiency comes with two 
drawbacks - the design of the correlator incurs large Non-Recurring Engineering (NRE) 



Frequency domain fractional sample correction could be performed by shortening the accumulation 
time, or the equivalent fractional sample correction could be performed in the time domain via a convolu- 
tion, but neither approach is generally computationally viable 



4-2. Correlator implementations 



51 



costs, and upgrading or altering the completed correlator is difficult (although this second 
point is mitigated somewhat in FPGA-based correlators). 

An alternative to a hardware correlator is to implement the correlation algorithm 
in a programming language such as C or C++, and run the correlation program in a 
generic multi-processor computing environment. Such a system will be hereafter referred 



to as a "software corre^ 



in software ( Bare et al. 



ator" . The first VLBI cor relators of the 1960's were implemented 



1967 



Moran et al 



19671 ) , but data processing requirements soon 



outgrew the capabilities of existing mainframe computers, and dedicated hardware systems 
were quickly developed. In recent years, the pace of advances in commodity computing 
hardware have begun to outstrip the demand increases of VLBI data processing, and 
software correlators have again become economically feasible. The first modern, high 
data-rate software correlator was developed for VLBI by the Com munications Research 
Laboratory (CRL) in Japan in the late 1990s (jKondo et al.l . |200J). Specialist software 
correlators have also been developed for purposes such as processing VLBI observations 
of spacec raft, which were use d to track the Huygens probe as it entered the atmosphere 



of Titan (Avruch et al 



20061). 



Compared to a hardware correlator, a software correlator offers several advantages: 
chiefly the ease and speed of development (minimising NRE), the adaptability of the sys- 
tem, and the flexibility of making alterations after the correlator has been deployed. The 
correlation algorithm is "embarrassingly parallel" and very well suited to commodity par- 
allel computing architectures, such as Beowulf clusters with Gigabit ethernet interconnect. 

Despite their obvious advantages in flexibility, by the mid-2000's no software corre- 
lators had yet been developed which were compatible with all the major VLBI formats 
currently used world-wide. As discussed in Chapter 1, at this time the LB A was complet- 
ing an upgrade which would provide the increased sensitivity necessary to undertake the 
pulsar astrometry program planned for this thesis, but which required a new correlator. 
Accordingly, a completely new, general-purpose software correlator was developed as part 
of this thesis. The primary goal of this correlator was to enable the LB A to undertake 
high-sensitivity observations and improve operational flexibility, but a secondary aim was 
for the correlator to be usable for VLBI (and indeed connected element interferometry) 
by many groups worldwide. The resultant "DiFX" software correlator is described in Sec- 
tion 14.31 befow, and its testing and verification are described in Section 14. 4[ Performance 
benchmarking of the DiFX code is described in Section [4.51 and some of the additional ap- 
plications it has enabled, beyond the pulsar astrometry of this thesis project, are described 
in Section I 
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4.3 DiFX 

DiFX (Distributed FX) is the name given to a code project encompassing a FX-style 
correlator (mpifxcorr) implemented in C++ and utilising message passing conforming to 
the Message Passing Implementation (Mpjfl) standard, as well as related functionalities 
such as geometric model generation (gencalc_delays), a graphical front-end (DiFXGUI) 
and various data and visibility inspection tools and plug-in packages. The DiFX package 
can be downloaded from http://astronomy.swin.edu.au/~adeller/software/difx/. 

DiFX requires an operating environment with a functional installation of MPI, and 
an implementation of a vector library. It was originally developed on an Intel cluster 
and by default makes use of the Intel Performance Primitives (IPlfl) library for vector 
calculations. The IPP library is is optimised for modern 32 and 64 bit CPU architectures 
and allows code acceleration of up to an order of magnitude compared to unoptimised 
code. It can run on any Intel or AMD CPU, but will only dispatch optimised code on 
Intel CPUs. An equivalent AMD library, the AMD Performance Library (APlfl) IS now 
available, but has not been tested. 



4.3.1 The DiFX code 

Figure 14.21 shows the high-level class structure of DiFX, along with the data flow. The 
correlation is managed by a master node (FxManager), which instructs data management 
nodes (Datastream) to send time ranges of baseband data to processing nodes (Core). 
The data are then processed by the Core nodes, and the results sent back to the FxMan- 
ager. Double buffered, non-blocking communication is used to avoid latency delays and 
maximise throughtput. Both the Datastream and Core classes can be (and have been) ex- 
tended to allow maximum code re-use when handling different data formats and processing 
algorithms. The Core nodes make use of an allocatable number of threads to maximise 
performance on a heterogenous cluster. 

The Datastream nodes can read the baseband data into their memory buffers from 
a local disk, a network disk or a network socket. Once the data are loaded into the 
datastream buffer, the remainder of the system is unaware of its origin. This is one of 
the most powerful aspects of this correlator architecture, meaning the same correlator can 
easily be used for disk-based VLBI correlation and real-time eVLBI, where the data is 
transmitted in real time from the telescopes to the correlator over optical fibre. The use 

*http://www-unix. mcs.anl.gov/mpi/ 

^http: / /www. intel.com /cd / software / products/asmo-na/eng/perflib/ipp / index. htm 
^http://developer. amd.com/apLlielp/aa_000_frames. html 
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Figure 4.2 Overview of the software correlator architecture. Data is loaded into memory 
from a disk or network connection by Datastream nodes. These nodes are directed by a 
Master node to send data from given time ranges (typically several ms) to the processing 
elements (Core nodes). The processed data are sent to the master node for long-term 
accumulation and storage on disk. 

of DiFX for eVLBI is discussed further in Section 14.6.11 

The implementation of the antenna- and baseline-based operations described in Sec- 
tion 14.11 that are necessary to form the visibility output of DiFX is described below. The 
additional processing required in the case of pulsar observations is also discussed. 

4.3.2 Antenna-based operations 
Alignment of telescope data streams 

To correlate data from a number of different telescopes, the changing delays between 
those telescopes must be calculated and used to align the recorded data streams at a 
predetermined point in space (in this case the geocentre) throughout the experiment. 

DiFX uses CALC ^ to generate a geometric delay model Tg{t) for each telescope in a 
given observation, at regular intervals (usually 1 second). CALC models many geometric 
effects, including precession, nutation, ocean and atmospheric loading, and is used by 
many VLBI correlators including the VLB A and EVN correlators. These delays are then 
interpolated (using a quadratic approximation) to produce accurate delays (Ar < 1 x 10^^^ 
seconds, compared to an exact CALC value) in double precision for any time during the 
course of the observation. The estimated station clock offets and rates are added to the 
CALC-generated geometric delays. 

^ http : // gemini .gsfc .nasa.gov /solve 
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The baseband data for each telescope are loaded into large buffers in memory, and the 
interpolated delay model is used to calculate the accurate delay between each telescope 
and the centre of the Earth at any given time during the experiment. This delay, rounded 
to the nearest sample, is the integer sample delay. The difference between the delay and 
the integer sample delay is recorded as the antenna based fractional sample delay (up to 
=b 0.5 sample). Note that the alignment of any two data streams (as opposed to a data 
stream alignment with the geocentre) is accurate to it 1 sample. 

The integer-sample delay is used to offset the data pointer in memory and select 
the data to be correlated (some number of samples which is an integer number of FFT 
windows, starting from the time of alignment). The fractional sample error is retained 
to correct the phase as a function of frequency following alignment to within one sample, 
fringe rotation, and channelisation (Section I4.3.2|) . 

Once the baseband data for each telescope have been selected, they are transferred 
to a processing node ("Core") using a non-blocking MPI send, and unpacked from the 
coarsely quantised representation (usually a 2-bit representation) to a floating point (sin- 
gle precision) representation. Presently supported VLBI formats includ e the Australia n 



standard LB ADR (Phillips et al., in preparation), Mark5A and Mark5B (jWhitnevl . 120031 ) , 



and the Japanese format K5 (jKovama et al.l . |2004| ) . If supported by the baseband format. 



missing or corrupted data is detected during the unpacking, and a count of valid samples 
is maintained. Ancillary data such as the sampled geometric delays at the beginning of 
each FFT window are also transferred (in double precision). From this point on, all oper- 
ations in the correlator are performed using floating point arithmetic, in single precision 
unless otherwise specified. Note that the data volume is expanded by a factor of 16 at this 
point. The choice of single precision floats (roughly double the precision necessary) was 
dictated by the capabilities of modern CPUs, which process single-precision floating point 
numbers efficiently. Using sufficient precision also avoids the small decorrelation losses 
incurred by optimised, low precision operations often used in hardware correlators. This 
is a good example of the sacrifice of efficiency for simplicity and accuracy with a software 
correlator. 

At this point all data streams from all telescopes are aligned to within it 1 sample 
of each other and the fractional sample errors for each of the telescope data streams are 
recorded for later use. A set number of samples from each telescope data stream have 
been selected and are awaiting unpacking and processing on a "Core" node (e.g. a PC in 
a Beowulf cluster). 
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Fringe rotation 

As shown in Section 14.1.31 for "pre-F" fringe rotation, the necessary complex fringe ro- 
tation function is assembled for each time sample by taking the sine and cosine of the 
geocentric delay multiplied by the sky frequency uq in a vector operation; it is applied via 
a vector complex multiplication for each telescope's data stream. The geometric delay at 
each time sample is obtained by interpolation between values supplied for the start and 
end of each FFT segment. 

Since the baseband data have already been unpacked to a floating point representation 
by this stage, a floating point fringe rotation is applied which yields no fringe rotation 
losses, compared, for example, to a 6.25% loss of signal to no ise for three level digital 
fringe rotation in a two level complex correlator ( Robertsl . 1997 ). 

Implemented as such, fringe rotation represents a mixing operation and will result in 
a phase difference term which is quasi-stationary at zero phase (the desired term) and a 
phase sum term which has a phase rate of twice the fringe rotation function, ~ 47rz/oT(t). 
The sum term vector averages to a (normally) negligible contribution to the correlator; 
for typical VLBI fringe rates (100s of kHz) and integration times (seconds) the relative 
magnitude of the unwanted contribution to each visibility point is < 10~^. In a software 
correlator it would be simple to control the integration time so that the rapidly varying 
phase term is integrated over exactly an integral number of terms of phase, thus making 
no contribution to the correlator output. This feature is not currently implemented in 
DiFX. 

As noted in Section [4.1.31 alternative method of fringe rotation for sufficiently low 
fringe rates is to apply a single correction per FFT window ("post-F" fringe rotation). In 
this case, no modulation is applied to the time domain data, but a single fringe-rotation 
value (appropriate for the midpoint of the FFT window) is calculated and retained for 
application during fractional-sampled correction. Post-F fringe rotation is desirable in 
situations where the fringe rate is extremely low, when the double-frequency term intro- 
duced by the mixing operation of pre-F fringe rotation is not effectively averaged to zero 
over the course of an integration and makes a significant and undesirable contribution to 
the correlator output. Switching from pre-F to post-F fringe rotation would be beneficial 
for periods of time in most experiments when the source traverses periods of low phase 
rate. Sources near a celestial pole can have very low fringe rates for long periods of time. 
Alternatively, if very short correlator integration times are used, the sum term may not 
integrate to zero when using pre-F fringe rotation. Post-F fringe rotation would therefore 
be a natural choice in these circumstances. 
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It should be noted that it is possible to undertake the exact equivalent to pre~F fringe 
rotation in the frequency domain. However, this would involve the Fourier transform of 
the fringe rotation function and a convolution in the frequency domain, which is at least 
as computationally intensive as the complex multiplication of the data and fringe rotation 
in the time domain. 

DiFX implements pre-F or post-F fringe rotation as a user controlled option. 

Channelisation and fractional sample error correction 

Once the data are aligned (and phase corrected, if pre-F fringe rotation is used), the 
time series data are converted into frequency series data (channelised), prior to cross 
multiplication. 

Channelisation of the data can be accomplished using an FFT. If pre-F fringe rotation 
has been applied, the data are already in complex form, and so a complex-to-complex 
FFT is used. The positive or negative frequencies are selected in the case of upper or 
lower sideband data respectively. If post-F fringe rotation is to be applied, the data are 
still real and so a more efficient real-to-complex FFT may be used. This is possible due 
to the conjugate symmetry property of an FFT of a real data series. In this case, lower 
sideband data may be recovered by reversing and conjugating the resultant channels. 

The final station-based operation is fractional-sample correction. The fractional sam- 
ple error is assumed not to vary over the FFT window, which is equivalent to the assump- 
tion made for post-F fringe rotation, but is considerably less stringent since the phase 
change is proportional to the bandwidth, rather than sky frequency as in the case of fringe 
rotation. As shown in Section [4.1.31 the frequency domain correction manifests itself as a 
slope in the phase as a function of frequency across the observed bandwidth. 

Thus, after channelisation, a further vector complex multiplication is applied to the 
channels, correcting the fractional sample error. In the case of post-F fringe rotation, the 
saved fringe rotation value is added to the fractional-sample correction for each frequency 
channel and the two steps are performed together. 

4.3.3 Baseline-based operations 

Cross multiplication of telescope data streams 

For each selected baseline, the channelised, compensated data from the telescope pair 
are cross-multiplied on a channel by channel basis (after forming the complex conjugate 
for the channelised data from one telescope), using a complex vector muliplication, to 
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yield the frequency domain complex visibilities that are the fundamental observables of 
an interferometer. This is repeated for each selected band/polarisation on each selected 
baseline. If dual polarisations have been recorded for any given band, the cross-polarisation 
terms can also be multiplied, recovering polarisation information for the target source. 

Integration of correlated output 

Once the above cycle of operations has been completed, it is repeated and the resulting 
visibilities accumulated (complex added) until a set accumulation time has been reached. 
Generally, on each cycle the input time increment is equal to the corresponding FFT length 
(twice the number of spectral points), but it is also possible to overlap FFTs. This allows 
more measurements of higher lags and greater sensitivity to spectral line observations, at 
the cost of increased computation. At the Core node, this "Short Term Accumulation" 
(STA) is equal to the number of FFT windows sent by the Datastream nodes. The 
STA results are then returned to the FxManager via a non-blocking MPI send, along 
with ancillary information such as the number of "good" samples for each baseline. The 
desired number of STA results are then accumulated at the FxManager node to reach the 
desired integration time. 

Calibration for nominal telescope Tgys 

Cross multiplication, accumulation and normalisation gives the complex cross power spec- 
trum for each baseline, representing the correlated fraction of the geometric mean of the 
powers detected at each telescope. To obtain the correlated power in units of Jy, the cross 
power spectra (amplitude components) should be scaled by the geometric mean of the 
powers received at each telescope measured in Jy i.e. the Tsys in Jy routinely measured 
at each antenna. Since the autocorrelation spectra are calculated concurrently with the 
cross-correlation spectra, this correction can be made within the correlator, before writing 
the visibilities to disk. Calibration is applied with a vector multiplication to each array of 
visibility spectra. Alternatively, the cross-correlation spectra can be simply be normalised 
by the number of contributing samples, and the calibration of the data can be completed 
offline, once measured, rather than a priori, Tgyg information is available. The application 
of normalisation or a priori system calibration is a user selectable option in DiFX. 

When the Tsys for each telescope is applied to calibrate the visibilities, it is also nec- 
essary to apply a scaling factor to compensate for decorrelation due to the coarse quan- 
tisation. This corrects the visibility amplitudes, but of course cannot recover the lost 
signal to noise. For the 2-bit data typically processed, this scaling factor is 1/0.88 in the 
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low-correlation limit (jCooperl . 119701 ) . The relationship becomes non- linear at high cor- 
relation and the scaling factor approaches unity as the correlation coefficient approaches 
unity. When applied on-line at the correlator, the low-correlation regime is assumed, and 
a correction for high-correlation cases can be applied in post-processing if necessary. 

Export of visibility data 

Once an accumulation interval has been reached, the visibilities must be stored in a useful 
format. DiFX can produce RPFIT^I format data, or an intermediate binary format which 



can be translated offline into F 



TS-IDn RPFITS and FX' 



'S-IDI ffles can be loaded into 



for data reduction. Ancillary 



analysis packages such as Aipf^l, CASA0, or MIRIAeEh 
information is included in the FITS ffle along with the complex visibilities, time stamps, 
and {u,v,w) coordinates. Presently, only FITS-IDI output supports the storage of visi- 
bility weights, which are adjusted based on the count of "good" samples detected in the 
unpacking process. 

4.3.4 Special processing operations: pulsar binning 

As discussed in Chapter^ pulsed signals are dispersed as they travel through the interstel- 
lar medium (ISM), resulting in a smearing of the pulse arrival time in frequency. In order 
to correct fo r the dispersive effects of the ISM, DiFX employs incoherent dedispersion 



( Vouteetal 



2002). This allows the visibilities generated by the correlator to be divided 
into pulse phase bins. Unlike hardware correlators which typically allow only a single 
on/off bin, or else employ 2^ bins of fixed width, DiFX allows an arbitrary number of bins 
placed at arbitary phase intervals. The individual bins can be written out separately in 
the FITS file format to enable investigation of pulse phase dependent effects, or can be 
filtered within the correlator to maximise signal to noise, based on a priori pulse profile 
information. 

To calculate which phase bin a visibility at a given frequency and time corresponds to, 
the software correlator requires information on the pulsar's ephemeris, which is supplied 
in the form of one or more "polyco" files containing a polynomial description of apparent 
pulse phase as a function of time. These are generated using the pulsar analysis program 
TEMPO0, and require prior timing of a pulsar. Additional software has been produced 



*http: //www. atnf.csiro.au/computmg/software/rpfits. html 
^http: //www. aoc.nrao.edu/aips/FlTS-lDI.fitml 

^"http: //www. aoc.nrao.edu/aips 

^^http://casa.nrao.edu/ 

^■^http: //www. atnf.csiro.au/computing/software/miriad 
^^http: / /pulsar. princeton.edu / tempo /reference_manual.litml 
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to verify the pulsar timing by using the generated polyco files and the baseband data from 
an experiment, allowing phase bins to be accurately set before correlation. 

For VLBI observations of pulsars, it is usually desirable to maximise the signal to noise 
of the observations by binning the visibilities based on the pulse phase, and applying a 
filter to the binned output based on the signal strength in that phase. Typically this filter 
is implemented as a binary on/off for each phase bin. Using the pulse profile generated 
from the baseband data of an observation, however, DiFX allows a user-specified number 
of bins to be generated and a filter applied based on pulse strength x bin width, allowing 
the maximum theoretical retrieval of signal, as described below. This also reduces the 
output data volume, since only an "integrated on-pulse" visibility is retained, rather than 
potentially many phase bins. 

Consider observing a single pulse, divided into M equally spaced phase bins. Let the 
pulsar signal strength as a function of phase bin be S{m), and the noise in single phase 
bin he Z X \/M, where Z is the baseline sensitivity for an integration time of a single pulse 
period. When all bins are summed (effectively no binning), the S/N ratio will be: 

Em=o S{m) 
Z 

as the signal adds coherently while the noise adds in quadrature. For a simple on/off gate 
accepting only bins mi to m2, the S/N ratio will be: 

^("^^ . (4.11) 



(Z X 



Finally, for the case where each bin is weighted by the pulse signal strength in that bin, 
the S/N ratio will be: 

E™=o(5M)' 



E^J^=o{Sim)xZx^]' 



(4.12) 



For a Gaussian shaped pulse, this allows a modest improvement in recovered signal to 
noise of 6% compared to an optimally placed single on/off bin. On a more complicated 
profile, such as a Gaussian main pulse with a Gaussian interpulse at half the amplitude, 
the improvement in recovered signal to noise increases to 21%. Throughout the remainder 
of this thesis, the term gate, when referring to a DiFX-correlated pulsar experiment, refers 
to a matched pulse filter in the sense described above. 
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4.3.5 Operating DiFX 

DiFX is controlled via an interactive Graphical User Interface (GUI), which calls the 
various component programs and helper scripts. The primary purpose of the GUI is to 
facilitate easy editing of the text files which configure the correlator, run external programs 
such as the delay model generator, and provide feedback while a job is running. Two files 
are necessary to run the actual correlator program. The first is an experiment configuration 
file, containing tables of stations, frequency setups, etc, analogous to a typical hardware 
correlator job configuration script. The second file contains the list of compute nodes on 
which the correlator program will run. 

While it is possible to run all tasks required to operate the correlator manually, in prac- 
tise they are organised via the GUI. This consists of running a series of helper applications 
from the GUI to generate the necessary input for the correlator. These include a script 
to extract experiment information from the VLBI exchange (VEX) file used to configure 
and schedule the telescopes at observe time, a delay and {u,v,w) generator which makes 
use of CALC 9, and scripts to extract the current load of available nodes. Pulsar-specific 
information such as pulse profiles and bin settings can also be loaded. This information 
is presented via the GUI and adjustments to the configuration, such as selecting compu- 
tational resources to be used, can be made before launching a correlation job. 

In the future it is planned to incorporate some real-time feedback of amplitude, phase 
and lag information from the current correlation via the GUI. This would be similar to the 
visibility spectra displays available continuously at connected-element interferometers. 

Alternative interfaces for configuring and launching DiFX have been developed by other 
users of the software, but are generally site-specific and interface to existing correlator 
control software. The foremost example of this approach is the VLB A, which is described 
in Section [4.6.21 



4.4 Deployment and verification 



In order to verify that the DiFX correlator was functioning correcting, three separate 
comparisons were made against existing, well-verified hardware correlators. In each case, 
the comparison showed complete agreement once differences in correlator architecture were 
accounted for. These comparisons are detailed below. 
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4.4.1 Comparison to the LB A S2 correlator 

Observations to provide data for a correlator comparison between DiFX and the ATNF 
S2 correlator were undertaken on March 12, 2006, with the following subset of the LBA: 
Parkes (64 m), ATCA (phased array of 5 x 22 m), Mopra (22 m), Hobart (26 m). 

Data from these observations were recorded simultaneously to S2 tapes and the LBADR 
disks (Phillips et al. 2008, in preparation) during a 20 minute period, UT 02:30-02:50, 
corresponding to a scan on a bright quasar (PKS 0208—512). Two 16 MHz right circular 
polarisation (RCP) bands were recorded, in the frequency ranges 2252 — 2268 MHz and 
2268 - 2284 MHz. 



The data recorded on S2 tapes were shipped to the ATNF LBA S2 correlator (j Wilson et al 

lOod ) at ATNF headquarters and processed. The S2 correlator is an XF-style hardware 



correlator. The data recorded to LBADR disks were shipped to the Swinburne University 
of Technology supercomputer and processed using the software correlator. 

At both correlators identical Tgyg values in Jy were specified for each antenna and 
applied in order to produce nominally calibrated visibility amplitudes, and both correlators 
used identical clock models, in the form of a single clock offset and linear rate as a function 
of time per antenna. The data were processed at each correlator using a two second 
integration time and 32 spectral channels across each 16 MHz band. 

Different implementations of the CALC-based delay generation were used at each cor- 
relator, and thus small differences exist in the applied delay models, which leads to differ- 
ences in the correlated visibility phase. These delay model differences were calculated and 
the phase due to differential delay model has been subtracted in the following discussion. 

From both cor relators, RPFITS format data were output and loaded into the MIRIAD 
software package (jSault et al.l . Il995l ) for inspection and analysis. The data from the two 
correlators are compared in Figures 14.31 - H31 

Figure 14.31 shows the visibility amplitudes for all baselines from both correlators as a 
function of time, over the period 02:36:00 - 02:45:00 UT, for one of the 16 MHz bands (2252 
— 2268 MHz). These amplitudes represent the vector averaged data over the frequency 
channel range 10 — 21 (to avoid the edges of the band). The data for each baseline were 
fit to a first order polynomial model {S{t) = ^t + 5o, where S is the flux density in Jy, t 
is the offset in seconds from UT 02:40:30, and Sq is the extrapolated flux density at time 
UT 02:40:30, using a standard linear least squares routine. The root mean square (RMS) 
variation around the best fit model was calculated for each baseline. The fitted models 
are shown in Figure |3]3] and show no significant differences between the S2 correlator and 
the software correlator. Further, the calculated RMS for each baseline agrees very well 
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Figure 4.3 S2 (red) and DiFX (black) visibility amplitude vs time for the 2252 - 2268 
MHz band on the source PKS 0208-512, as described in the text (PKS = Parkes; MOP = 
Mopra; HOB = Hobart; NAR = ATCA). Symbols represent the actual visibilities produced 
by the correlators, while the lines represent linear least-squares fits to the visibilities (one 
line per dataset). 



between DiFX and the S2 correlator, as summarised in Table IT2l 

Figure 14.41 shows the visibility phase as a function of time for each of the six baselines 
in the array. Again the data represent the vector averaged correlator output over the 
frequency channel range 10 — 21 within the 2252 — 2268 MHz band. As discussed above, 
small differences between the delay models used at each correlator have been taken into 
account as part of this comparison. 

Figure 14.51 shows a comparison of the visibility amplitudes and phases as a function of 
frequency in the 2252 — 2268 MHz band. The data represented here result from a vector 
average of the two datasets over a two minute time range, UT 02:40:00 — 02:42:00. Since 
the S2 correlator is an XF-style correlator, it cannot exactly correct fractional sample error 
in the same manner as an FX correlator such as DiFX, as the channelisation is performed 
after accumulation. The coarse (post-accumulation) fractional sample correction leads to 
decorrelation at all points except the band center, up to a maximum of ~ 10% at the band 
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Figure 4.4 S2 (red) and DiFX (black) visibility phase vs time for the 2252 - 2268 MHz 
band on the source PKS 0208—512, as described in the text. Antenna labels are identical 
to Figure 133] above. The PKS-NAR baseline has been shifted by —50° for clarity. 



Table 4.2. Linear fit parameters for visibility amplitude vs time for DiFX and the LBA 

S2 correlator, with 95% confidence limits 



Baseline 


OffsetDiPX (Jy) 


OffsetLBA (Jy) 


SlopeDiFX (pJy s ^) 


SlopCLBA (/Jjy s ^) 


PKS - NAR 


1.341 ± 0.030 


1.343 ± 0.028 


10 ± 13 


14 ± 12 


PKS - MOP 


3.185 ± 0.058 


3.185 ± 0.063 


14 ± 24 


-11 ± 26 


PKS - HOB 


2.307 ± 0.058 


2.293 ± 0.061 


-12 ± 24 


- 6 ± 24 


NAR - MOP 


1.616 ± 0.109 


1.619 ± 0.114 


-27 ± 43 


-10 ± 45 


NAR - HOB 


1.142 ± 0.111 


1.139 ± 0.116 


- 3 ± 44 


- 5 ± 46 


MOP - HOB 


2.694 ± 0.256 


2.681 ± 0.257 


18 ± 101 


56 ± 101 
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Figure 4.5 S2 (red) and DiFX (black) visibility amplitude and phase vs frequency data 
for the 2252 - 2268 MHz band on the source PKS 0208-512, as described in the text. 
Antenna labels are identical to Figure 14.31 above. The S2 data has been corrected for 
fractional-sample error decorrelation at the band edges as described in the text. 
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edges on long baselines where the geometric delay changes by a sample or more over an 
integration period. Compensation has been performed for this band edge decorrelation in 
the S2 correlator amplitudes shown in Figure 



4.4.2 Comparison to the VLB A correlator 

Data obtained as part of a regular series of VLBA test observations were used as a basis 
for a correlator comp arison between DiF X and the VLBA correlator, which is an FX-style 
hardware correlator ([Napier et al.l . Il994l ). The observations were made on August 5, 2006 
using the Brewster, Los Alamos, Mauna Kea, Owens Valley, Pie Town, and Saint Croix 
VLBA stations. One bit digitised data sampled at the Nyquist rate for four d ual polari- 
satio n bands, each of 8 MHz bandwidth, were recorded using the Mk5 system (jWhitneyl . 



20031 '). The four bands were at centre frequencies of 2279.49, 2287.49, 2295.49, and 2303.49 



MHz. The experiment code for the observations was MT628 and the source observed was 
0923+392, a strong and compact active galactic nucleus. Approximately two minutes of 
data recorded in this way was used for the comparison. 

The Mk5 data were correlated on the VLBA correlator and exported to FITS format 
files. The data were also shipped to the Swinburne supercomputer and correlated using 
DiFX, producing RPFITS format files. In both cases, no scaling of the correlated visibility 
amplitudes by the system temperatures were made at the correlators. The visibilities 
remained in the form of correlation coefficients for the purposes of the comparison - i.e. 
a system temperature of unity was used to scale the amplitudes. Each 8 MHz band was 
correlated with 64 spectral points, and an integration time of 2.048 seconds was used. 

The VLBA correlator data were read into AIPS using FITLD with the parameter 
DIGIC0R=1. The DIGICOR parameter is used to apply certain scalings to the visibility 
amplitudes for data from the VLBA correlator. Further, to obtain the most accurate 
scaling of the visibility amplitudes, the task ACCOR was used to correct for imperfect 
sampler thresholds, deriving corrections to the antenna-based amplitudes of ~ 0.5%. These 
ACCOR corrections were applied to the data, which were then written to disk in FITS 
format. 

The software correlator data were read directly into AIPS and then written to disk 
in the same FITS format as the VLBA correlator data. No corrections to amplitude or 
phase of the software correlated data were made in AIPS. 

The VLBA correlator data and the software correlator data were both imported into 
MIRIAD for inspection and analysis, using the same software as used for the comparison 
with the LBA correlator described above. RCP from the 2283.49 - 2291.49 MHz band 
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Figure 4.6 VLBA correlator (red) and DiFX (black) visibility amplitude vs time for the 
2283.49 " 2291.49 RCP band from the VLBA test observation MT628, as described in the 
text. The units of time are seconds from UT 00:00:00, and the amplitude scale is correlation 
coefficient. Symbols represent the actual visibilities produced by the correlators, while the 
lines represent linear least-squares fits to the visibilities. The text annotation on each 
panel lists the average correlation coefficient amplitude for each correlator over the time 
period, as tabulated in Table 



over the time range UT 17:49:00 — 17:51:00 was used in all comparison plots below. 

Since the delay models used by the VLBA and software correlators differ at the pi- 
cosecond level, as is the case for the comparison with the LBA data in Section 14.4. H 
differences in the visibility phase exist between the correlated datasets. As with the LBA 
comparison, compensation has been performed for the phase errors due to the delay model 
differences in the following comparison. 

Figure 14.61 shows the visibility amplitudes for all baselines from both correlators as a 
function of time. These amplitudes represent the vector averaged data over the frequency 
channel range 10 — 55 (to avoid the edges of the band). The data for each baseline were fit 
to a first order polynomial model {S{t) = + Sq, where S is the correlation coefficient, t 
is the offset in seconds from UT 17:50:00, and Sq is the extrapolated correlation coefficient 
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Figure 4.7 VLBA correlator (red) and DiFX (black) visibility phase vs time for the 2283.49 
- 2291.49 RCP band from the VLBA test observation MT628, as described in the text. 
The units of time are seconds from UT 00:00:00, and phase is displayed in degrees. 



Table 4.3. Linear fit parameters for visibility amplitude (in units of correlation 
coefficient) vs time for DiFX and the VLBA correlator, with 95% confidence limits 
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Figure 4.8 VLBA correlator (red) and DiFX (black) visibility amplitude and phase as a 
function of frequency for the 2283.49 - 2291.49 RCP band from the VLBA test observation 
MT628, as described in the text. The vertical scale for correlation coefficient amplitude 
on each panel is - 0.018, while the phase scale spans ±180°. The horizontal scale for 
each panel displays channels 0-64. 



4.4- Deployment and verification 



69 



at time UT 17:50:00) using a standard least squares regression. The RMS variation around 
the best fit model was calculated for each baseline. The fitted models are shown in Figure 
14.61 and show no significant differences between the VLBA correlator and DiFX. Further, 
the calculated RMS for each baseline agrees very well between the VLBA correlator and 
DiFX. The results of the comparison are summarised in Table 14. 3[ 

Figure [121 shows the visibility phase as a function of time for each of the fifteen baselines 
in the array. Again, the data represents the vector averaged correlator output over the 
frequency channel range 10 — 55 within the band. As discussed above, small differences 
between the delay models used at each correlator cause phase offsets between the two 
correlators, and have been taken into account as part of this comparison. 

Figure 14.81 shows a comparison of the visibility amplitudes and phases as a function 
of frequency in the band. The data represented here results from a vector average of 
the two datasets over a two minute time range. Figures 14. 6t 14.71 and 14.81 show that the 
results obtained by the VLBA correlator and DiFX agree to within the RMS errors of the 
visibilities in each expected. Since the VLBA correlator is an FX-style correlator, 

like DiFX, there is no requirement to correct for band-edge decorrelation as was required 
for the LBA correlator. 



4.4.3 Comparison to the MPIfR geodetic correlator 

Unlike the S2 and VLBA correlators, the MarklV correlator operated by the Max Planck 
Institut fiir Radioastronomie (MPIfR) in Bonn, Germany, is used prim arily for geodeti c 
observations. The MarklV correlator, an XF-type hardware correlator (jWhitnevl . Il993l ). 
has been used heavily for geodetic VLBI processing for many years as part of the Interna- 
tional VLBI Service (IVS; http:/ /ivscc. gsfc.nasa.gov/) and is one of a very small number 
of trusted geodetic correlators around the world. 

The comparison with this correlator is described in detail in lTingay et al.l (j2008l ). and 
the results are summarised here. The data used for the correlator comparison was a subset 
of the data collected for a geodetic experiment (BM261) conducted with the NRAO VLBA 
on 2007 July 03. One minute of data (11:23:00 UT - 11:24:00 UT) from four antennas 
were selected for the purposes of the comparison: Fort Davis (FD), Pie Town (PT), 
Owens Valley (OV), and Kitt Peak (KP). 8 dual polarisation bands of bandwidth 8 MHz 
were recorded, with the following centre frequencies: 1350.55, 1358.55, 1366.55, 1374.55, 
1382.55, 1390.55, 1398.55, 1406.55 MHz. For the comparison, the RCP data from the 
1350.55 MHz band was used. 

The data at the VLBA antennas were recorded to Mark5 disk packs and transported 
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to the MarklV correlator at MPIfR, where correlation was performed with the MarklV 
correlator. The comparison subset of data was then exported from the MarkS units to a 
standard filesystem and correlated using DiFX on a commodity cluster of linux machines 
at MPIfR. The data were correlated with 128 frequency channels across the band. Prior 
to correlation, a clock model was derived for the experiment and applied identically for 
each correlator. As with the previous two comparisons, however, geometric models were 
derived independently, and the model difference was calculated and subtracted from the 
following analysis. 

Figure 14.91 shows the visibility phase as a function of time for each of the six baselines 
correlated. As with the previous comparisons, the data represents the vector averaged 
correlator output from the centre of the band (frequency channels 32 — 96) , and the phase 
offsets due to small differences between the delay models used at each correlator have been 
corrected in this analysis. Figure H. 101 shows a comparison of the visibility amplitudes and 
phases as a function of frequency in the band. The data represented here results from a 
vector average of the two datasets over the one minute comparison time range. 

Figures 14.91 and 14.101 show that the results obtained by the MarklV correlator and 
DiFX are in good agreement. Since the MarklV correlator is used primarily for geodesy 
and the amplitude correction signal path is not well documented, no attempt has been 
made to improve the a priori amplitude corrections in the manner of the VLBA correlator 
comparison. As the MarklV correlator is an XF-style correlator and the visibilities were 
not corrected for fractional-sample decorrelation at the band edge, a noticeable difference 
can be seen at the band edge in amplitudes. Figure 14.91 also shows several glitches in 
visibility phase in the MarklV correlator output. The origin of these glitches is unknown, 
but they are almost certainly an error in the MarklV correlator output, as no glitches 
are seen to occur on correlators simultaneously and the magnitude of the phase jump is 
many times the phase RMS. A possible cause of the phase jumps is an error on the MarkS 
Station Unit during playback at the correlator (C. Phillips, private communication). 

4.5 Performance benchmarking 

In order to keep every compute node used in the correlation fully loaded, they must 
be kept supplied with raw data. If this condition is satisfied, the correlation is CPU- 
limited, and the addition of further nodes will result in a linear performance gain. In 
practise, however, at some point obtaining data from the data source (network socket or 
disk) and transmitting it across the local network to the processing nodes will no longer 
occur quickly enough, and the correlation becomes data-limited rather than CPU-limited. 
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Figure 4.9 Geodetic correlator (red) and DiFX (black) visibility amplitude vs time for 
the 1346.55 - 1354.55 RCP band from the geodetic observation BM261, as described in 
the text. The units of time are seconds from UT 00:00:00, and the phase is displayed in 
degrees. 
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Figure 4.10 Geodetic correlator (red) and DiFX (black) visibility amplitude and phase as 
a function of frequency for the 1346.55 - 1354.55 RCP band from the geodetic observation 
BM261, as described in the text. The vertical scale for correlation coefficient amplitude 
on each panel is - 0.05, while the phase scale spans ±180°. The horizontal scale for each 
panel displays channels 0-128. 
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Correct selection of correlation parameters, and good cluster design, will minimise the 
networking overhead imposed on a correlation job, and ensure that all compute nodes are 
fully utilised. This is discussed in Section 14.5.11 below, and performance profiles for the 
CPU-limited case are presented in Section 14.5.21 

4.5.1 Networking considerations 

As described in Section 14.3.11 double-buffered communications to the processing nodes 
are used to ensure that nodes are never idle as long as sufficient aggregate networking 
capability is available. The use of MPI communications adds a small but unavoidable 
overhead to data transfer, meaning the maximum throughput of the system is slightly less 
than the maximum network capacity on the most heavily loaded data path. 

There are two significant data flows: out of each Datastream and into the FxManager. 
For any high speed correlation, there will be more Core nodes than Datastream nodes, so 
the aggregate rate into a Core will be lower than that out of a Datastream. The flow out 
of a Core is a factor of A'^corcs times lower than that into the FxManager node. 

If processing in real time (when processing time equals observation time), the rate 
out of each Datastream will be equal to the recording rate, which can be up to 1 Gbps 
with modern VLBI arrays and is within the capabilities of modern commodity ethernet 
equipment. The rate into the FxManager node will be equal to the product of the recording 
rate, the compression ratio, and the number of Cores, where the compression ratio is the 
ratio of data into a Core to data out of a Core. This is determined by the number of 
antennas (since number of baselines scales with number of antennas squared), the number 
of channels in the output cross-power spectrum, the number of polarisation products 
correlated, and the integration time used before sending data back to the FxManager 
node. 

It is clearly desirable to maximise the size of data messages sent to a core for pro- 
cessing, since this minimises the data rate into the FxManager node for a given number 
of Cores. However, if the messages are too large, performance will suffer as RAM ca- 
pacity is exceeded. Network latency may also become problematic, even with buffering. 
Furthermore, it should be apparent that in this architecture, the Cores act as short-term 
accumulators (STAs), with the FxManager performing the long term accumulation. The 
length of the STA sets the minimum integration time. It is important to note, however, 
that the STA interval is entirely conflgurable in the software correlator, to be as short as 
a single FFT, although network bandwidth and latency are likely to be limiting factors in 
this case. 
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For the majority of experiments it is possible to set a STA length which satisfies all the 
network criteria and allows the Cores to be maximally utilised. For combinations of large 
numbers of antennas and very high spectral and time resolution, however, it is impossible 
to set an STA which allows a satisfactorily low return data rate to the FxManager node. 
In this case, real time processing of the experiment is not possible without the installation 
of additional network and/or CPU capacity on the FxManager node. 

It is important to emphasise that although it is possible to find experimental configu- 
rations for which DiFX suffers a reduction in performance, these configurations would be 
impossible on existing hardware correlators. If communication to the FxManager node is 
limiting performance, it is also possible to parallelize a disk-based experiment by dividing 
an experiment into several time ranges and processing these time ranges simultaneously, 
allowing an aggregate processing rate which equals real time. This is actually one of the 
most powerful aspects of the software correlator, and one which would allow scheduling of 
correlation to always ensure the cluster was being fully utilised. 

4.5.2 CPU-limited performance 

Figure 14.111 shows the results of performance testing on the Swinburne cluster available in 
2006 (3.2 GHz Pentium 4 machines, Gigabit ethernet network) for different array sizes and 
spectral resolutions. The results shown in Figure H.lll were obtained for data for which the 
aggregate bandwidth was 64 MHz, broken up into 8 bands each of 8 MHz bandwidth (4 x 
dual polarisation 8 MHz bands: data were 2-bit sampled: antenna data rate 256 Mbps). 
Node requirements for real-time operation are extrapolated from the compute time on an 
8 node cluster. The correlation integration time is 1 second and all correlations provide 
all four polarisation products. RAM requirements per node ranged from 10 - 50 MB 
depending on spectral resolution, showing that large amounts of RAM are unnecessary 
for typical correlations. It can be seen that even a modestly sized commodity cluster can 
process a VLBI-sized array in real time at currently available data rates. 

Recently, the Swinburne supercomputer cluster has been upgraded to multicore, 64 bit 
machines. Testing on this architecture has shown that the correlation code performance 
scales near-linearly with the number of threads used per node (up to the number of 
CPU cores), as expected for this embarrassingly parallel problem. Data rates of up to 1 
Gbps per antenna have been correlated in real-time on these new machines, in the eVLBI 
experiments described in Section [4.6.11 



4-6. Additional applications 



75 



Aggregate bandwidth: 64 MHz 
Number of channels: 256 
Integration time: 1 s 



10 



15 



20 



Number of antennas 



o 















Aggregate bandwidth: 


64 MHz 






Number of antennas: 


10 : 






ntegratlon time: 1 s 






■Si- 


<■ 





10 100 1000 10^ 10^ 



Number of channels 

Figure 4.11 Benchmark data showing the computational requirements of DiFX to corre- 
late in real-time, as described in the text. The nodes are single core 3.2 GHz Pentium 
processors with 1 GB RAM, and in both benchmarks 64 MHz of total bandwidth per 
station was correlated with a 1 second integration period. Top panel shows the scaling 
of computational requirements with number of antenna, using 256 spectral points per 8 
MHz subband. Bottom panel shows the scaling of computional requirements with spectral 
points per subbband for a ten station array. 



4.6 Additional applications 

Whilst DiFX was developed with the primary aim of servicing the needs of the LBA, every 
effort was made during development to ensure that DiFX would be adaptable to the needs 
of other interferometers that required an upgrade path from an existing correlator. This 
approach has enabled the use of DiFX to correlate data from every major VLBI array in 
the world, as well as facilitating the rapid development of several new instruments. Those 
science highlights which are external to this thesis are summarised briefly below. 

4.6.1 LBA eVLBI 

Operational eVLBI modes have been tested using DiFX, using real-time data from the 
three ATNF telescopes (Parkes, ATCA, and Mopra), the University of Tasmania telescope 
in Hobart, and international telescopes including the Kashima 34m telescope in Japan, 
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and the Shanghai 25m telescope in China. These observations have seen data rates up 
to 1 Gbps per antenna using Austrahan-only antennas, and 512 Mbps per antenna in 
international experiments, equaling the highest data rate experini ents performed with the 
EVN (highest published data rate 256 Mbpili; iParagi et al.l . 120071 ) . Correlation has taken 



place using computing resources at Swinburne University, the ATNF, and the University 
of Western Australia in Perth (a Cray XD-1 utilising Opteron processors and on-board 
Xilinx FPGAs). A m ore extensive descrip tion of the eVLBI activities of the LB A using 



DiFX can be found in 



PhiUips et alJ (|2007|). 



First science results from these eVLBI experiments have rec ently been published , and 
include observations of the Galactic x-ray bi r iary C ircinus X-1 (jPhillips et all 120071 ) and 
the supernova remnant of SN1987A (|Tingavl . 120081 ) . 



4.6.2 VLBA sensitivity upgrade 

The VLBA is the only full-time VLBI instrument in the world, and has been responsible for 
much of the rapid advan ce in VLB! scien ce over the past decade. When designed, it used 
the Mark4 tape svs tem (IWhitnevi. Il993l '). and it has been upgraded to use the Mark5A 



disk-based system ( Whitneyl . l2003l ). Recently, a major upgrade to VLBA sensitivity 
has been planned, aiming to reach data rates of 4 Gbps, an 8-fold increase over the 
current maximum data rate. Information on the planned sensitivity upgrade is available 
at http://www.vlba.iirao.edu/memos/sensi/. In order to reach these increased data 
rates, the VLBA requires both a new digital backend and a new correlator. The digital 
backend, named Mark5C, is under development, and DiFX has been chosen as the upgrade 
path for the correlator. "First light" with VLBA data on a local installation of DiFX 
has been obtained in 2008, and production usage is expected to commence in late 2008. 
Support and configuration tools have been developed to allow DiFX to function as a 
"drop-in" replacement for the existing hardware correlator, allowing duplication of existing 
functionality at the same time as providing new capabilities. 



4.6.3 High Sensitivity Array observations of pulsar scintillation 

As discussed in Chapter [21 the inhomogeneous distribution of the ISM leads to diffractive 
and refractive scattering of signals from compact sources such as pulsars. Whilst single- 
dish studies of pulsar scintillation have already shed considerable li ght on the nature of 
the structures in the ISM responsible for this scintillation (see e.g. 



Cordesetal 



20061 . 



*for unpublished higher data rate results, see http://www.expres-eu.org/512Mbps_6tel.html 
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and references therein), high frequency resolution VLBI observations of scintillating pul- 
sars offer the possibility of directly imaging the ISM structures at extremely high spatial 
resolution. 

The capabilities of DiFX for this type of analysis have recently been demonstrated 
with observations of PSR B0834— 04 in 2005 (Brisken et al., in preparation). Elements of 
the High Sensitivity Array (HSA), comprising the NRAO Green Bank Telescope (GBT; 
100 m), Westerbork (14 x 25 m), Jodrell Bank (76 m), and Arecibo (305 m) were used to 
provide an ultra-sensitive array at 327 MHz. The data were recorded using the Mk5 system 
and correlated on a prototype implementation of DiFX. The main requirement on the 
correlation was 250 Hz wide frequency channels, over the broadest bandwidth available, to 
maximise signal to noise. For these observations a 32 MHz band was available. The DiFX 
correlation used 131,072 frequency channels, obtaining a frequency resolution of 244 Hz. 
Such extreme frequency resolution for a continuum experiment is beyond the capabilities 
of any existing hardware correlator. Figure B. 121 shows a section of the dynamic spectrum 
from this observation which shows the scintillation structure as a function of time and 
frequency on the GBT-Arecibo baseline. A deconvolution process allows the generation of 
a speckle image of the scattering disk using the interfeometric phase information (Brisken 
et al., in preparation). The effective resolution of the speckle image is better than 1 
mas, revealing structure which is a factor of 40 smaller than the synthesised beam of the 
interferometer. 

4.6.4 EVN pulsar astrometry 

DiFX has been used to correlate Mark5 data from the EVN with multiple pulsar gates 
and phase centres to make an image of the pulsar distribution in the globular cluster M15 
(Deller et al., in preparation). Through phase referencing to a nearby quasar and taking 
an ensemble average of pulsar motions over several epochs, it will be possible to obtain 
the proper motion of the cluster with a higher signal/noise ratio than would be possible 
using a single pulsar. It is also possible to further improve the calibration of the data 
by using the brightest pulsar as an in-beam calibrator, improving the sensitivity to the 
other pulsars in the cluster and making it possible to conduct extremely precise relative 
astrometry to probe the intra-cluster dynamics. 

4.6.5 Widefield VLBI 

DiFX is being used to correlate VLBA observations of the Chandra Deep Field South, using 
high time and frequency resolution which is not available with the existing VLBA hardware 



78 



Chapter 4. DIFX: AN FX STYLE SOFTWARE CORRELATOR 




Figure 4.12 The cross-power dynamic spectrum showing scintihation variations for the 
pulsar B0834— 04 on the Green Bank Telescope - Arecibo baseline. Brightness represents 
the visibility amplitude and colour represents the visibility phase. Increasing frequency 
runs left to right and increasing time runs top to bottom. This section of the dynamic 
spectrum represents just 5% of the time span and 0.5% of the bandwidth of the observation 
(330 seconds and 160 kHz). 

correlator. This project aims to obtain VLBI-resolution imaging over a 30 x 30 arcminute 
field of view. Post-correlation phase shifting to potential targets identified from lower- 
resolution observations allows many targets to be imaged from a single correlator pass. In 
the future, DiFX will support multiple phase centres during correlation, eliminating the 
need for post-correlation shifting and averaging of visibilities before imaging, which will 
enable this and other larger-scale VLBI surveys to be carried out quickly and efficiently. 



4.6.6 Australian/New Zealand geodesy 

A new 3-element, full time geodetic array is planned to commence operations in Australia 
in 2009 as part of the AuScope projeclo. This array will contribute to observations which 
address both local and global geodetic goals. Closely linked to the AuScope project is the 
development of a geodetic VLBI capacity in New Zealand, led by the Centre for Radio- 
physics and Space Science at Auckland University of Technologjo. The DiFX software 
correlator will be used for the correlation of these new geodetic arrays. 



"'^^http://www. auscope.org.au/ 

^^http://www.aut.ac.nz/about/faculties/design_and_creative_technologies/faculty_research_ofBce/ 
centre_for_radiophysics_and_space_research/mission.htm 



4-6. Additional applications 79 
4.6.7 Worldwide geodesy 

As is the case with many areas of VLBI, the geodetic community is currently planning a 
significant upgrade in capabilities. This upgrade is spearheaded by a sensitivity increase 
driven by wider recording bandwidths, but also includes new hardware such as smaller, 
more rapidly slewing antennas and new, improved digital backends. The upgrade project 
is known as VLBI2O1C0, and aims to achieve an initial recording rate of 2-4 Gbps with 
the potential for expansion to 8, 16 or even 32 Gbps. As the Mark 4 hardware correlators 
currently used for geodetic VLBI are limited to a 1 Gbps correlation rate, th is project will 
necessitate new and improved correlator infrastructure ( Niell et al. . 2OO5I ). At present. 



the use of DiFX in this role is being investigated by MPIfR, with activities including the 
correlator comparison shown in Section [4.4.31 



http://ivs.nict.go.jp/mirror/about/wg/wg3/ 



5 

ASTROMETRIC OBSERVATIONS, DATA 
REDUCTION AND ANALYSIS 



5.1 Goals and target selection 



The observational program of this thesis encompassed 8 pulsars, which are listed in Ta- 
ble 15.11 As already noted in Chapter [Tl previous Southern Hemispher e VLBI pulsar as- 



trometry programs (jPodson et al 



200 



I 



Leggd, 



2002 



Bailes et al 



1990) have yielded only 



two published pulsar parallaxes, so one primary motivation was to increase the number of 
southern VLBI parallaxes and hence improve models of Galactic electron distributions at 
southern declinations. Additionally, as this is the first such large scale southern parallax 
study, target pulsars of varying brightness and predicted distance were chosen in order to 
determine the types of targets that would be feasible for future southern VLBI studies. 
The brightness and predicted distance of the selected targets are shown below in Table I^TTl 
Within the bounds of these criteria, however, it was possible to choose target pulsars for 
which a VLBI parallax would enable a deeper understanding of additional astrophysical 
phenomena. 

The chosen targets can be divided into three broad categories: binary pulsars used for 
tests of GR and gravitational wave detection (PSR J0437-4715, PSR J0737-3039A/B and 
PSR J2145-0750), pulsars with unusual luminosity in the radio (low luminosity pulsars 
PSR J0108-1431 and PSR J2144-3933) or x-ray (PSR J0630-2834, whose x-ray luminosity 
is anomalously high), and "technique check" sources which are bright and predicted to 
be nearby (PSR J1559-4438 and PSR J2048-1616). The results for each pulsar, and the 
implications of these results, are discussed individually in Chapter [6l 

The remainder of this chapter will be devoted to describing the observations and data 
reduction applied to all pulsars, with a particular focus on the tools developed for pulsar 
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Table 5.1. Target pulsars 



Pulsar 
name 


DAI distance 
(pc)" 


1600 MHz 
flux (mjy) 


Pulsar 
gating gain 


Equivalent gated 
1600 MHz flux (mJy) 


Reference 
source 


Calibrator / target 
separation (deg) 


J0108-1431 


130 


0.6 


5.1 


3 


JOlll-1317 


1.5 


J0437-4715 


140 


140 


6.25 


875 


J0439-4522 


1.9 


J0630-2834 


2150 


23 


3.5 


81 


J0628-2805 


0.7 


J0737-3039A/B 


570 


1.6 


2.5 


4 


J0738-3025 


0.4 


J1559-4438 


1600 


40 


3.6 


144 


J1604-4441 


0.9 


J2048-1616 


640 


13 


3.6 


47 


J2047-1639 


0.5 


J2144-3933 


180 


0.8'' 


10 


8 


J2141-3729 


2.1 


J2145-0750 


500 


8" 


4.3 


34 


J2142-0437 


3.3 



^Takcii from Taylor &i Cordcs (1993) 

^Pulsar suffers heavily from long timescale scintillation, so individual epochs vary considerably from the average value shown. 

astrometry reduction which address aspects of the data reduction process not previously 
highlighted in the literature. In order to illustrate the concepts, the results for PSR 
J1559-4438 will be used, as this proved the superior of the two "technique check" sources. 

5.2 Observations 

Eight observational epochs were spread over a two year period from May 2006 to February 
2008. Epochs were typically 24 hours in duration and a subset of the eight pulsars were 
observed, depending on which were closest to parallax extrema. The observing frequency 
was centred on 1400 MHz for the first observation, and 1650 MHz for the remaining seven 
observations. PSR J0437-4715 was observed separately, with four epochs of 12 hours 
duration, centred on 8400 MHz. Observations at this higher frequency were made possible 
by the high flux and narrow pulse profile of PSR J0437-4715 . Dual circular polarization 
was used at all frequencies. 

The Australian Long Baseline Array (LBA) consists of six antennas - the Australia 
Telescope National Facility (ATNF) telescopes in New South Wales (Parkes, Australia 
Telescope Compact Array [ATCA], Mopra); the University of Tasmania telescopes at 
Hobart, Tasmania and Ceduna, South Australia; and the NASA DSN facility at Tid- 
binbilla, Australian Capital Territory. The Parkes, phased ATCA, Mopra, and Hobart 
telescopes participated in all experiments, but a Tidbinbilla antenna (70m or 34m) was 
used only when available, and Ceduna participated in observations of J0437-4715 onl}0. 
The maximum baseline length with Ceduna is 1700 km, and without Ceduna is 1400 km. 
Representative uv coverage at 1650 MHz and 8400 MHz is shown in Figure EHJ 

All observations used the recently introduced LBADR disk-based recording system 

'^Ceduna does not possess a 1600 MHz receiver, so could only participate in the higher frequency 
experiments 
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5 -5 4D 20 -20 -40 

U (10'*) " ('0°« 

Figure 5.1 Typical uv coverage at 1650 MHz (no Ceduna; left panel) and at 8400 MHz 
(with Ceduna; right). The target sources are PSR J1559-4438 and PSR J0437-4715 
respectively. Without Ceduna, the uv coverage is heavily biased north-south, and wide 
hour-angle coverage is necessary to gain acceptable uv coverage. The displayed plots are 
optimal uv coverage - equipment failure or telescope commitments often meant only a 
subset of this coverage was obtained. 



(Phillips et al., in preparation). At the three ATNF observatories, the presence of two 
Data Acquisition System (DAS) units allowed a recording rate of 512 Mbps (8 x 16 
MHz bands, Nyquist sampled at 2 bits), while the non-ATNF stations recorded at 256 
Mbps (4 X 16 MHz bands). For epochs where dual-polarization feeds were available 
at all antennas, two frequency bands were dropped at the non-ATNF stations, but as 
some of the NASA DSN feeds are single polarization only, 1650 MHz epochs featuring 
the 70m NASA antenna and 8400 MHz epochs featuring the 34m NASA antenna instead 
retained a single polarization of all frequency bands. The observation dates, targets, and 
participating antennas are summarised in Table | 



A phase reference cycle time of six minutes, apportioned equally to target and calibra- 
tor, was used for all observations. As the LBA consists of disparate antennas ranging up to 
70m in diameter (and, when phased, the ATCA has an equivalent diameter of hundreds of 
metres for the purposes of calculating field of view), it was not possible to utilise in-beam 
calibrators for ariy sources, u nlike recent pulsar astrometric programs using the VLBA 



(jChatteriee et al.l . I2OOII . l2005|) . 



The data were correlated using matched filtering on pulse profiles with the DiFX soft- 
ware correlator, producing RPFITS format visibility data. Table ISTT] shows the predicted 
gain due to pulse profile filtering for each target source. Pulsar ephemerides were ob- 
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Table 5.2. Observation summary 



Observation Duration Observed pulsars Participating telescopes 

date (MJD) (hours) 



53868 


12 






0437-4715 




Parkcs, ATCA, Mopra, Hobart, Ccduna 


53870 


15 


1559- 


-4438, 2048-1616, 2144-3933, 2145- 


-0750 


Parkos, ATCA, Mopra, Hobart 


53970 


24 




0108- 
1559- 


-1431, 0630-2834, 0737-3039 
-4438, 2048-1616, 2144-3933 




Parkos, ATCA, Mopra, Hobart, DSS43 (70m) 


54055 


12 






0437-4715 




Parkes, ATCA, Mopra, Hobart, Cedunaf" DSS43 (70m 


54057 


24 




0108- 


-1431, 0630-2834, 0737-3039 




Parkes, ATCA, Mopra, Hobart, DSS43 (70m) 






1559- 


-4438, 


2048-1616, 2144-3933, 2145- 


-0750 




54127 


24 




0737- 


0108-1431, 0630-2834 
-3039, 1559-4438, 2144-3933 




Parkos, ATCA, Mopra, Hobart, DSS43 (70m) 


54181 


12 






0437-4715 




Parkes, ATCA, Mopra, Hobart, Ceduna, DSS43 (70m) 


54182 


24 




0630- 
2048- 


-2834, 0737-3039, 1559-4438 
-1616, 2144-3933, 2145-0750 




Parkes, ATCA, MoprapHobart, DSS43 (70m) 


54307 


24 




0737- 


0108-1431, 0630-2834 
-3039, 1559-4438, 2144-3933 




Parkes, ATCA, Mopra, Hobart, DSS43 (70m) 


54413 


24 




0108- 


-1431, 0630-2834, 0737-3039 




Parkes, ATCA, Mopra, Hobart 






1559- 


-4438, 


2048-1616, 2144-3933, 2145- 


-0750 




54417 


12 






0437-4715 




Parkes, ATCA, Mopra, Hobart, Ceduna, DSS34 (34m) 


54500 


24 




0737- 


0108-1431, 0630-2834 
-3039, 1559-4438, 2144-3933 




Parkes, ATCA, Mopra, Hobart, DSS43 (70m)° 



Equipment failure meant no useful data was available from Ccduna for this session. 
^Limited time {~ 4 hours) was available at the DSS telescopes in these sessions. 
^Equipment failure led to the loss of half of the Mopra data from this session. 
^Setup errors meant no useful DSS43 data was available from this session. 



tained using the ATNF Pulsar Catalogue (jManchester et al.l . |2005| ) , with the exception of 
the double pulsar PSR J0737-3039A/B, which was periodically updated with the latest 
published ephemeris. Two second integrations and 64 spectral points per 16 MHz band 
were used for all observations. 



5.3 Data reduction 



Data re duction was pe r forme d principally in AIPS, usin g the pyt h on in terface Parsel- 



Tongue (jKettenis et al 



The DIFMAP package (jShepherdl . 119971 ) was used for 
imaging and self calibration. The data reduction was implemented as a pipeline, with 
user interaction for imaging, editing of solution tables, and visibility flagging. All scripts 
used in the data reduction process are freely available and can be downloaded from 
http : / / astronomy . swin . edu . au/~adeller/ software/scripts/. The individual stages 
of the pipeline are described below. 



5.3.1 Amplitude and weight calibration and flagging 

Amplitude calibration using the measured telescope system temperatures was carried out 
using the AIPS tasks APCAL and ANTAB. Flagging based on predicted or (when avail- 
able) logged telescope off-source times due to slewing or failures was applied using UVFLG, 
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while the first and last 10 seconds of every scan was excised with the task QUACK. For 
the ATCA, the tied array infrastructure required flagging the first three correlator inte- 
grations (totaling 30 seconds) of each scan with QUACK. The weight of each visibility 
point, which is set to a constant value when the RPFITS format data is loaded into AIPS, 
was initially scaled by the predicted baseline sensitivity using a ParselTongue script. The 
effect of data weighting is investigated further in Section 15.61 



5.3.2 Geometric model and ionospheric corrections 

At the low frequencies which are generally used for pulsar astrometry, ion ospheric varia- 
tions usually make the dominant contribution to systematic error (see e.g. 



Brisken et al, 



2OO2I ) . Using ionospheric models based on Global Position System (GPS) data provided 



by the NASA Jet Propulsion LaboratorjB, the AIPS task TECOR was used to correct 
phase variations due to the ionosphere. This observing program roughly coincided with 
the solar minimum of 2006, and consequently the ionospheric variations were generally 
at a minimum. For PSR J0437-4715, which was observed at 8.4 GHz, the position shifts 
introduced by ionospheric correction were very small (~100 /uas) and did not improve the 
quality of the astrometric fit, so this step was not applied for PSR J0437-4715. 

Total Electron Content (TEC) models suffer at southern declinations due to the rela- 
tively sparse distribution of GPS receivers at southern latitudes. Consequently, the derived 
ionospheric corrections are much less reliable for the LBA than for similar Northern Hemi- 
sphere instruments. The dispersive delay corrections generated by TECOR were inspected 
for each epoch, and the effectiveness of different TEC maps is investigated in Section [5.5.21 

While the observational program was underway, considerably more accurate station 
positions were derived for several LBA an tennas using archival geodetic observations and 
the OCCAM software ( Titov et al. . l2004 V a dedicated 22 GHz LBA geodetic experiment 



(Petrov et al., in preparation), and GPS survey measurements. Additionally, more ac- 
curate positions for some ca librators were published in the 5th VLBA Calibrator Survey 



(VCS5; iKovalev et al.l . 120071 ). The visibilities were corrected to account for the revised 



positions using an AIPS SN table generated with the Wizardry feature of ParselTongue, 
which stored the difference between the initial and corrected geometric models. 

For some pulsars, their proper motions (> 100 mas yr~^) cause significant position 
shifts over the course of a 24 hour observation, comparable in some cases to the epoch 
positional accuracy. As the geometric model generation used in DiFX at the time of these 



^available from the Crustal Dynamics Data Information Systems (CDDIS) archive: 
ftp://cddis.gsfc.nasa.gov/pub/gps/products/ionex/ 



86 



Chapter 5. OBSERVATIONS, DATA REDUCTION AND ANALYSIS 



observations could not account for proper motion, the visibility phases and uvw values 
were corrected in AIPS using a SN table generated by ParselTongue, interpolating between 
predicted postions for the pulsar at the start and end of the experiment. This also allowed 
the proper motion to be refined (or, in some cases, measured for the first time) over the 
course of the observations before the final corrections were applied. 

5.3.3 Fringe-fitting and amplitude calibration refinement 

Fringe fitting was performed using the AIPS task FRING, using a point source model, on 
the phase reference calibrator data for each target pulsar. Subsequently, a single structural 
model was produced for each calibrator using the combined datasets from all epochs. With 
the exception of B0736~303, each source was modeled using a dominant component (delta 
or narrow Gaussian) fixed at the phase center, and 0-2 secondary components which 
were allowed to vary in positiorlfl. The flux density of all components was allowed to vary. 
Typically, the variation in flux density of secondary component (s) between epochs was 
< 1% of the image peak flux density. Images were then generated using uniform weighting 
(bin size two pizels) and baseline weights set to predicted noise RMS (difmap uvweight 
settings 2,-1). The images of each phase reference source are shown in Figures 15.21 and 
15.31 With the exception of B0736-303 and J2142-0437, the corrections due to reference 
source structure were very small. 

The solutions were applied to each calibrator and the data averaged in frequency, 
exported to disk and loaded into DIFMAP. The calibrator model was loaded and several 
iterations of self-calibration and modelfitting performed. The difmap 'modelfit' command 
uses the Levenberg-Marquardt least-squares minimisation algorithm to fit the free model 
parameters to the visibility points, incorporating the visibility weights. The self calibration 
corrections were then written to disk as an AIPS SN table using the 'cordump' patch to 
difmapQ. Additionally, data points flagged in DIFMAP were collated in a Wizardry script 
and converted into a flag file suitable for the AIPS task UVFLG. 



^B0736-303 was modeled using a combination of CLEAN and modelfit components, and is discussed 
further in Section [6.1.21 

*http://astronomy.swin. edu.au/~elenc/difmap-patclies/ 
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Figure 5.2 LBA images of the phase reference sources JOlll-1317, J0439-4522, J0628- 
2805 and B0736-303. 0.4 mas pixels were used for J0439~4522 (8 GHz observations) and 
2 mas pixels used elsewhere. Uniform weighting with visibility weights raised to the power 
— 1 was used for all images. 
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Figure 5.3 LBA images of the phase reference sources J1604-4441, J2047-1639, J2141- 
3729 and J2142-0437. 2 mas pixels and uniform weighting with visibility weights raised 
to the power —1 were used for all images. 
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The amplitude corrections generated in this manner ahowed compensation for imper- 
fectly measured system temperature values, and were also applied to visibility weights. 
For bands which could not be self-calibrated, due to only being observed by the three 
ATNF antennas, a correction was estimated based on the nearest band with self calibra- 
tion solutions of the same polarization. The self calibration solutions were loaded into 
AIPS using TBIN, and applied to the target pulsars using CLCAL. 

The use of bandpass calibration was investigated but found to make insignificant dif- 
ference to the fitted target position. The LBA Data Acquisition System (DAS) utilitizes 
digital filtering and typical bandpass phase ripple was < 2 degrees. When averaging in 
frequency, the lowest and highest 10% of the band was discarded and the central 80% of 
the band averaged with uniform weight assigned to each channel. 



5.3.4 Pulsar scintillation correction 



Nearby pulsars can suffer dramatic, and rapid, variations in visibility amplitude due 
to diffractive scintillation. The size of the scintillation pattern is typically much larger 
than the size of the Earth, and so the amplitude variations are essentially independent of 
baseline length. Maximal amplitude fluctuations (where the RMS is equ al to the mea n 
flux density) are seen for pulsars in the strong scattering regime (see e.g. IWalkeri . Il998l ). 
which can be pred i cted f rom Galactic electron distribution models. The NE2001 model 



( Cordes and Lazio 



2OO2I ) predicts that strong scintillation should be observed for PSR 



J0630-2834, PSR J1559-4438, PSR J2048-1616, PSR J2144-3933 and PSR J2145-0750. 
An example is shown for PSR J1559-4438 in Figure 15.4b . which shows the variation of 
visibility amplitude with time for Tidbinbilla baselines over a 2 hour period. Observed 
scintillation parameters for target pulsars are shown in Table 15. 3[ Assuming the material 
responsible for the scintillation is turbulent, with a Kolmogorov distribution, the scintilla- 
tion time Tscint and scintillation bandwidth -Bgdnt can be scaled to t he frequencies used in 
these observations with the relations Tsdnt oc u^'"^ and -Bscint oc z/^-^ ( Cordes et al. . 19861 ). 

Left uncorrected, the variation in visibility amplitude with time scatters power through- 
out the image plane, as shown in Figure [5^ . which shows the residuals for PSR J1559- 
4438 after fltting a single point source to the visibilities shown in Figure Eilt^- To overcome 
this, a ParselTongue script was written to produce an AIPS SN table which would flatten 
visibility amplitudes by averaging data over all sensitive baselines to 1 /4 of the scintillation 
time, normalizing, and taking the square root to obtain an antenna based correction. The 
visibility weights were then scaled by the inverse of the square of the correction, upweight- 
ing points when the amplitude was high and downweighting points of low signiflcance. The 



90 Chapter 5. OBSERVATIONS, DATA REDUCTION AND ANALYSIS 




RIgnl AilcEnsien fmuB^J Right Ascension fmas) 



Figure 5.4 The effects of diffractive scintillation for PSR J1559~4438. Panel a) shows the 
uncorrected visibility amplitude on baselines including the Tidbinbilla antenna from one 
experiment - the scintillation timescale of several minutes is apparent. Panel b) shows 
the same visibility amplitudes after correcting for scintillation. Panel c) shows the image 
residuals of the uncorrected dataset - contours are 1,2,4 and 8 mJy/beam. Panel d) shows 
image residuals after correcting for scintillation, with contours at 1,2 and 4 mJy/beam - 
the improvement in image quality is obvious. 



5.4- Positional determination and parallax fitting 
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Table 5.3. Observed pulsar scintillation parameters and estimated scattering disk sizes 



Pulsar 


Scintillation 


Scintillation 


Observing 


Reference 




Scattering disk size^ 


name 


time (s) 


bandwidth (MHz) 


frequency (MHz) 






(mas, 1C50 MHz) 


J0108-1431 




* 


* 


Johnston et al. 


(1998) 


0.003 


J0437-4715 


600 


17.4 


660 


Johnston et al. 


(1998) 


0.052 


J0630-2834 


456 


0.2 


327 


Bhat et al. 


(1999) 


0.022 


J0737-3039A 


> 70'^ 


1.4'' 


1400 


Coles et al. 


(2005) 


0.406 


J1559-4438 


77 


0.16 


660 


Johnston et al. 


(1998) 


0.133 


J2048-1616 


138 


0.54 


327 


Bhat et al. 


(1999) 


0.023 


J2144-3933 


1500 


2.9 


660 


Johnston et al. 


(1998) 


0.126 


J2145-0750 


1510 


1.48 


436 


Johnston et al. 


(1998) 


0.046 



^Calculated from dccorrclation bandwidth where available, and taken from NE2001 model where scintillation measurements 
are unavailable 

'^Scintillation parameters have not been measured, but are believed to be extremely large, as expected for a very nearby 
pulsar 

'^Varies considerable with orbital velocity — values quoted for highest transverse velocity 
*^Coles, private communication 

effect of these corrections is shown in Figures [5.4b and l5.4d . which show the visibility am- 
plitudes and image residuals respectively. In this example, the signal to noise ratio (SNR) 
of the detection is improved by a factor of two when the visibility amplitudes are corrected 
for scintillation, which reduces the position determination error by a corresponding factor. 

5.4 Positional determination and parallax fitting 

The calibrated visibilities for each pulsar were averaged in frequency, written to disk and 
loaded into DIFMAP. A single delta function model component was initialized at the 
peak of the dirty image, and the modelfit command used to obtain the best fit for the 
pulsar visibility data. Each observing band, as well as the combined dataset, was then 
imaged separately using uniform weighting (as before with a bin size of 2 pixels and 
weights set to predicted baseline SNR) producing 8 images which were saved and read 
back into AIPS using the task FITLD. Variations to the weighting scheme are discussed 
in Section 15.5.31 The AIPS task JMFIT was used to determine the pulsar position and 
formal (SNR-based) errors in the image plane. Systematic offsets of several (3-6) mas 
were observed at all epochs between the 4 bands that were only contributed to by the three 
ATNF antennas, and the 4 bands common to all antennas. The magnitude and direction 
of these offsets varied between epochs, and so the ATNF-only bands (which additionally 
had formal errors of approximately five times the other bands, primarily due to the shorter 
baselines) were dropped. It is suspected that these offsets were caused by the lack of an 
accurate amplitude self-calibration solution for these bands. 
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For each pulsar, the best fit values of J2000 position (right ascension and declination), 
proper motion (right ascension and declination) and parallax were initially determined 
by iteratively minimising the error function calculated by summing the value of predicted 
minus actual position over all measurements, weighted by t he individual measur ement 



errors. The iterative minimisation code used is described in iBrisken et al.l (j2002l ). and 
accounts for the covariances between parameters when calculating the uncertainty of fitted 
parameters. A reference time for the proper motion was chosen to be 31 Dec 2006 (MJD 
54100) to minimize proper motion uncertainty contibutions to be position uncertainty. 

This approach yields error estimates for the 5 fitted parameters which are almost 
certainly an underestimate, for two reasons: 

1. There are systematic errors, varying from band to band within an epoch (intra-epoch 
errors) such as the residual unmodeled differential ionosphere, bandpass effects etc, 
and varying from epoch to epoch (inter-epoch errors) caused by effects dependent on 
observing time, such as seasonal or d iurnal ionosph eric variations, and variations in 



refractive scintillation image wander (jRickett 



3). 



19901 ). These should increase the error 



on each individual measurement, but estimating the magnitude of the systematic 
component for each individual measurement is poorly constrained. 

2. Each measurement is assumed to be completely independent, whereas as noted 
above correlated errors are expected between measurements from the same epoch. 
In essence, this approach overestimates the number of independent measurements, 
lowering the reduced chi-squared and implying a better fit than the actual result. 

It is possible to make an estimate of the magnitude of the intra-epoch errors by comparing 
the scatter in fitted positions from the individual bands. Forming a weighted centroid 
position for each epoch utilising all measurements from that epoch allows an estimation 
of the likelihood that the measured points are consistent with that centroid, through the 
calculation of a reduced chi-squared value. If the reduced chi-squared value exceeds unity, 
the presence of unmodeled systematic errors can be inferred. 

Since there is no a priori knowledge of the systematic error distribution, an equal 
systematic error was allocated to to each measurement, and added in quadrature to the 
original measurement error. The errors in right ascension and declination are treated 
separately. This is necessarily an iterative procedure, since the weighted centroid will 
be altered by the addition of these systematic error estimates. In effect, this assumes 
a zero mean, gaussian distribution for the systematic errors. Although this is unlikely 
to be the true distribution, it is the most reasonable assumption available, and certainly 
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more correct than assuming no systematic errors at all. The intra-epoch systematic error 
estimate for an epoch is obtained when the reduced chi-squared reaches unity. 

Once a single position measurement and error has been calculated for each epoch, the fit 
to position, parallax and proper motion was re-calculated, and the reduced chi-squared 
of the fit inspected again. Since the degrees of freedom were reduced, the addition of 
systematic errors did not always result in a lowering of reduced chi-squared. The presence 
of inter-epoch systematic errors was inferred if the reduced chi-squared remained greater 
than unity. Again, without knowledge of the underlying distribution of these errors, it was 
decided to apportion an equal amount to each epoch, iterating until a reduced chi-squared 
of unity was obtained. As with intra-epoch errors, right ascension and declination were 
treated separately. 

Thus, the final astrometric dataset for each pulsar consisted of a single position mea- 
surement for each epoch, with a total error equal to the weighted sum of the individual 
band formal errors, added in quadrature to the estimates of intra-epoch and inter-epoch 
systematic errors. 

For each pulsar, the robustness of error estimation after the inclusion of estimated sys- 
tematic contributions was checked by implementing a bootstrap technique. Bootstrapping, 
which involves repeated trials on samples selected with replacement from the population 
of measured position points, is a statistical technique allowing the estir nation of paramete r 
errors without a complete knowledge of the underlying distribution (jPress et al.l . 12002 ). 
It differs from Monte Carlo analysis in that bootstrapping uses the original sample set as 
a population to draw from, whereas Monte Carlo techniques generate data samples based 
upon assumed parameter distributions. In this instance, the original single band position 
measurements were taken as the population, and N samples were drawn, where N was the 
original number of measurements. Each bootstrap consisted of 10,000 such trials, and the 
parameter errors estimated from the variance of the resultant distributions. A minimum of 
3 different epochs needed to be included to ensure a fit was possible - on the rare occasion 
that a trial did not satisfy this requirement, it was re-drawn. 

Thus, three sets of fitted parameters and errors were obtained for each pulsar - a 
"naive" result using the single-band positions (which generally possessed a reduced chi- 
squared greater than 1.0), a bootstrap result, and a more conservative estimate which 
attempts to account for the impact of systematic errors (the "inclusive" fit). In general, 
the estimated magnitude of errors on fitted parameters increased through these three 
different schemes. Typically, the ratio in the errors on the inclusive fit to those on the 
naive fit ranged from 0.95 to 1.90. The final error values obtained from the inclusive fit 
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are the most accurate estimation possible, and are believed to be inherently conservative. 
All quoted errors are la unless otherwise stated. 

Pulsar positions are given in the ICRF frame, and position errors are the formal astro- 
metric fit errors only and do not include the component due to reference source position 
uncertainty. With the exception of B0736-303, all reference sources have been previously 
observed in ICRF or VLBA Calibrator Series (VCS) observations, and possess position 
uncertainties ranging from 0.8 mas to 1.6 mas. B0736-303 was observed using the VLBA 
in phase-referencing mode, with ICRF source B0736-332 used as a calibrator. This rela- 
tive astrometry obtained the position of B0736-303 in the ICRF frame, but the precision 
was limited to tens of milliarcseconds, due to the large differential ionospheric contribution 
at the low elevation of the sources for the VLBA. 

This position error is unacceptably large, and accordingly the position of B0736-303 
was "reverse-engineered" from the astrometric fit to the position of PSR J0737-3039A/B. 
Since the position of PSR J0737-3039A/B is known to better than one mas in the Solar 
System reference frame, and the Solar System frame and ICRF are aligned at the several 
mas level, the VLBI position obtained from the astrometric fit could be compared to the 
timing position and a calibrator offset deduced. As with the VLBA phase referencing 
to B0736-332, a small position offset due to phase referencing errors can be expected, 
but as PSR J0737-3039A/B is much closer to B0736-303 (~20 arcminutes, as opposed 
to three degrees for B0736-332), and is observed at higher elevation using the LBA, 
these residual errors are much smaller. As described in Section 16.1.21 the position of 
B0736-332 was adjusted to ensure a match between the timing and VLBI positions for 
PSR J0737-3039A/B, resulting in a final positional accuracy for B0736-303 estimated 
at approximately five mas. In hindsight, phase-referencing of B0736-303 to B0736-332 
should also have been performed using the LBA, but insufficient time was available for 
such an observation after this position error was discovered, and a single LBA observation 
would be unlikely to improve on the estimated five mas accuracy eventually determined 
for B0736-303. 

5.5 Technique check: PSR J1559-4438 

5.5.1 Initial results 

Using the techniques described above, initial results were obtained for PSR J 1559-4438 
(shown in Table [53]). The motion of PSR J1559-4438 in right ascension and declination 
is shown in Figure ESI along with the fitted path according to the inclusive fit. The fit is 
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Motion of pulsar J1 559-4438 




-5-4-3-2-101234 

Relative RA (mas) from 15^59'^41 .526077^ 



Figure 5.5 Motion of the pulsar in right ascension and dechnation, with measured positions 
overlaid on the best fit. Sensitivity-weighted visibilities were used. The motion of the 
pulsar is positive in right ascension and declination. The first epoch (lower left) is clearly 
inconsistent. 

clearly unsatisfactory, since it predicts a negative parallax (though consistent with zero). 
The final column of Table 15.41 shows that systematic errors far exceed the nominal single- 
epoch positional accuracies, and inspection of Figure [531 shows that the first epoch (MJD 
53870) is markedly discrepant with the remaining epochs. 

As shown below, fine-tuning of the data reduction is required in order to obtain op- 
timal results from each pulsar, in particular with regard to the details of the ionospheric 
corrections and the visibility weighting schemes using in imaging. The steps taken for 
PSR J1559-4438 in are described below in Sections 15.5.21 and 15.5.31 

5.5.2 Ionospheric correction 

The obvious source of the large systematic errors present in the initial fit shown in Sec- 
tion 15.5.11 is the varying ionospheric correction between epochs. Accordingly, as a first 
check, the position fits for each epoch were recalculated after subtracting the differential 
ionospheric correction between PSR J1559-4438 and its phase reference - in effect, remov- 
ing the applied ionospheric correction and leaving the data uncorrected for ionospheric 
effects. This was implemented using a ParselTongue script which made a two-point inter- 
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Table 5.4. Initial results (sensitivity-weighted visibilities) for PSR J1559-4438 



Parameter 



Naive fit 



Bootstrap fit 



Inclusive fit 



Right ascension (J2000) 
Declination {J2000) 
fia (mas yr~ ^ ) 
Ma (mas yr~^) 
7r (mas) 

Nominal distance (pc) 
Nominal vt (km s~ ) 
Reduced chi— squared 
Mean epoch fit error (mas) 
Intra— epoch sys. error (mas) 
Inter— epoch sys. error (mas) 



15:59:41.526092 ± 0.000005 
-44:38:45.902028 ± 0.000017 



2.64 ± 0.07 
13.36 ± 0.02 
-0.066 ± 0.038 
15100 

-974 
8.8 



15:59:35.526093 ± 0.000031 
-44:38:45.902034 ± 0.000062 



2.65 ± 0.44 
13.36 ± 0.06 

-0.054 ± 0.193 

-18500 
-1190 



15:59:41.526077 ± 0.000026 
-44:38:45.901989 ± 0.000102 
3.08 ± 0.46 
13.29 ± 0.17 
-0.083 ± 0.245 



11900 
-770 
1.0 
0.198 
0.097 
1.070 



polation between adjacent calibrator scans to calculate the differential correction to the 
target (which was not absorbed into the fringe-fit), which was stored in a CL table and 
subtracted using the AIPS task SPLAT. The applied values were saved for later analysis. 

Intuitively, the largest ionospheric corrections would be expected when the angular 
displacement of the pulsar from the sun is small, since this angular separation (along with 
the level of solar activity) is the largest influence on the TEC. The angular displacement at 
each epoch between the original fitted position and the position obtained when ionospheric 
correction was removed is presented in Table [53| and plotted against angular separation of 
the pulsar from the sun at the time of observation in Figure 15.61 The revised astrometric 
fit obtained without ionospheric correction is plotted in Figure E31 

It is immediately apparent from Figure 15.61 that the first epoch (MJD 53870) is dis- 
crepant in that the position shift due to ionospheric correction is unusually large, given 
the large angular separation from the Sun. As shown in Figure 15.71 this epoch becomes 
more consistent with the fit when the ionospheric correction is removed, this epoch be- 
comes more consistent with the fit, but the third epoch (MJD 54057) becomes much more 
inconsistent. This is unsurprising, however, since this epoch had the smallest pulsar-Sun 
separation and the largest ionospheric corrections. 

To investigate whether the chosen ionospheric map was at fault, the first epoch was 
re-reduced with all available maps from the NASA CDDIS archivqj, but no significant 
change was found in fitted position. Given that the different TEC maps make use of many 
of the same GPS stations, this is unsurprising. This problem is exacerbated at southern 
declinations due to the low density of GPS receivers at southern latitudes. Additionally, 
any errors in the TEC maps would have an impact ~40% greater for this first epoch, 
due to its lower observing frequency of 1400 MHz, compared with the 1650 MHz center 
frequency used for all subsequent observations. 

^ftp: / /cddis. gsfc.nasa.gov / pub /gps / products /ionex / 
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Table 5.5. Average position shift due to ionospheric corrections for PSR J1559-4438 



Epoch MJD 


RA shift (mas) 


Dec shift (mas) 


Sun separation (deg) 


53870 


-5.36 


0.27 


153.5 


53970 


-1.34 


0.15 


96.4 


54057 


-6.61 


0.39 


25.8 


54127 


-3.15 


-0.07 


63.0 


54182 


-2.66 


0.35 


113.6 


54307 


-2.51 


0.34 


120.1 


54413 


-3.68 


0.20 


30.8 


54500 


-2.89 


0.05 


69.8 



• 1650 MHz epochs 
■ 1400 MHz epoch 



:0 40 60 80 100 120 140 160 

Target - sun angular separation (deg) 



Figure 5.6 Shift in fitted position due to ionospheric correction for PSR J1559~4438 vs 
angular separation between the pulsar and the Sun. The single 1400 MHz epoch has an 
unusually large correction given the large angular separation between the pulsar and the 
Sun. 
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Figure 5.7 Motion of the pulsar in right ascension and dechnation, with measured positions 
overlaid on the best fit, when ionospheric corrections have been removed. Sensitivity- 
weighted visibilities were used. The first epoch (lower left) is now more consistent, but 
the third epoch (during which the pulsar-Sun angular separation was only 26 degrees) is 
now inconsistent. 



Thus, due to the probability of residual ionospheric errors for this epoch, and also 
the potential for frequency-dependent calibrator source structure, the first epoch (MJD 
53870, the only 1400 MHz epoch) was dropped from all further analysis. The fit to the 
remaining seven epochs, with ionospheric corrections re-enabled, is shown in Figure 15.81 
While a realistic fit is now obtained, the measurement of parallax is still not significant. 
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Motion of pulsar J1 559-4438 




-3-2-1 1 2 3 



Relative RA (mas) from 15^59'^41.526138^ 

Figure 5.8 Motion of the pulsar in right ascension and dechnation, with measured positions 
overlaid on the best fit, with ionospheric corrections reinstated but the first epoch dropped. 
Sensitivity- weighted visibilities were used. The fit is improved considerably. 



100 



Chapter 5. OBSERVATIONS, DATA REDUCTION AND ANALYSIS 



5.5.3 Data weights 

Initial pulsar imaging and position fitting used visibilities weighted according to the best 
estimate of instantaneous baseline sensitivity. Whilst this is theoretically optimal for 
data which consists only of signal S and additive random thermal errors £'thcrmi it fails 
to account for the presence of multiplicative systematic errors E'sys = 6*'^''''=^*-' caused 
by residual calibration errors. Typically, these systematic errors are dominated by at- 
mospheric and ionospheric gradients, although other contributions include antenna and 
calibrator so urce positio r i erro rs, time-variation of antenna bandpasses, and instrumental 



phase 
while 



jitter. 



Foma 



Pradeletal 



ont 



20051 ) presents a theoretical review of phase referencing errors. 



(120061 ) presents a simulation-based approach. 
If (j)sys{t) had zero-mean and was ergodic, its effect would be indistinguishable from 
thermal noise and could be easily estimated, allowing the visibility weights to be cor- 
rected. For antenna/source position errors and large-scale atmospheric/ionospheric struc- 
ture, however, the residual errors are correlated over long times, causing systematic shifts 
in the fitted position for the target. 

When normal sensitivity-based weighting is employed in the presence of substantial 
and persistent systematic phase errors, the systematic noise on the most sensitive baseline 
will be absorbed into the fit, at the cost of a poorer fit to the less sensitive baselines. The 
magnitude of the induced error will be dependent on the ratio of systematic to thermal 
errors, and the discrepancy in sensitivity between the most and least sensitive baselines 
in the array. For the LBA, the most sensitive baseline (Parkes-DSS43: system equivalent 
flux density 30 Jy) is roughly 13 times more sensitive than the least sensitive (Hobart- 
Mopra: 380 Jy). Thus, the LBA is particularly susceptible to the influence of systematic 
errors, due to the pronounced variation in baseline sensitivities. 

The systematic errors can be crudely estimated (in a model-dependent fashion) by per- 
forming phase-only self-calibration on the target pulsar over a sufficiently long timescale 
to obtain sufficient SNR, and comparing the magnitude of the corrections to those expected 
from thermal noise alone. While this approach probes systematic errors over a somewhat 
shorter time period than those which would dominate for inter-epoch errors (tens of min- 
utes, rather than hours), it is illustrative of the presence of systematic errors overall. Such 
corrections are shown in Figure 15.91 for the ATCA station using a three-minute solution 
interval during the observation on MJD 54127. The corrections are clearly correlated over 
timescales of tens of minutes, and the RMS deviation of 4.4° is an order of magnitude 
greater than the estimated thermal phase RMS of 0.4° (calculated in this high-SNR limit 
as station sensitivity divided divided by target flux density, scaled by pulse filtering gain 
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Figure 5.9 Self-calibration corrections, using a three-minute timescale, for the ATCA 
station on PSR J1559-4438. Clear systematic deviations are seen from the zero-mean 
distribution expected from purely thermal noise. The RMS of the corrections exceeds 
those expected due to thermal noise by an order of magnitude. 



and converted from radians to degrees). Thus, for this observation, systematic errors ^ 
thermal errors and weighting visibilities by sensitivity actually degrades the quality of the 
position fit. 



If the average £^sys could be accurately estimated for each baseline over the duration 
of an experiment, the baseline visibility weights could be adjusted by assuming that Egys 
is time-invariant. Even more desirable would be the estimation of Egys as a function 
of time, allowing a time-variable adjustment of the visibility weights. Given present 
instrumentation, there is no way to reliably estimate systematic error (time dependent or 
independent) in a model independent fashion. In the limit where ii^sys ^ -E'thcmn however, 
the visibility weight for each baseline (regardless of sensitivity) will be dominated by the 
systematic error contribution, resulting in approximately equal weights for all visibilities. 
Accordingly, visibility weights for all baselines were reset to an equal, constant value and 
data reduction repeated. The results are discussed below in Section [5.5.4[ 
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5.5.4 Final results 

The revised fit obtained using equally weighted visibihty data is described in Table 15.61 
and plotted in Figure 15.101 Through comparison of Table 15.61 with Table 15.41 it can be 
seen that the average fit error for a single epoch has increased by 20%, but the inter-epoch 
systematic error has decreased by 95%. Thus, while using equally weighted data incurs 
a small sensitivity penalty, it benefits significantly through the reduced susceptibility to 
systematic errors. The fitted distance of 2600 pc is consistent with the NE2001 distance 
of 2350 pc, as well as HI absorption measurements which imply a lower distance limit of 



2000 lb 500 pc (jKoribalski et al.l . Il995l ). The implications of this distance measurement are 
discussed in Section 16.2.31 

The use of natural weighting, as opposed to uniform weighting, was investigated but 
found to produce inferior results. Fitted parameters remained relatively constant but 
errors on the fitted parameters increased by ~ 50%. This is unsurprising, since the use of 
natural weighting promotes a larger beamsize due to the concentration of more visibility 
points at small uv distances. Natural weighting is obviously more likely to prove useful 
when -Etherm ^ -^sys; and its use for the weakest target pulsars is shown in Chapter [6l 



5.6 Optimal data weighting 

From the results shown in Sections 15.5.31 and 15.5.41 it is clear that for PSR J1559-4438 the 
astrometric error budget is dominated by systematic errors, and that the use of equally 
weighted visibility points is optimal. However, Table [5TT] shows that this may not be the 
case for other pulsars in the target sample, as some are orders of magnitude fainter than 
PSR J1559-4438. Accordingly, the conditions under which sensitivity-weighted visibilities 
give superior results to equally-weighted visibilities were investigated. This was carried 
out by adding simulated thermal noise of varying RMS to the existing dataset. 

Three "noisy" datasets D^, Db, and Dc were constructed by adding zero-mean, 
gaussian-distributed noise to the real and imaginary visibility components of the original 
observations. Since the theoretical single epoch SNR for sensitivity- weighted data should 
be ~ 800 (a factor of 10 greater than the typically obtained SNR), the RMS of the added 
noise in the three datasets was set to predicted baseline sensitivity scaled by a factor of 20, 
40, and 80, which should allow a maximum single-epoch SNR of 40, 20 and 10 respectively. 
The results of fitting the modified datasets (using the inclusive fit approach only) with and 
without the use of sensitivity weighting are presented in Tables 15.71 and 15.81 respectively. 
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Figure 5.10 Motion of the pulsar in right ascension and declination, with measured posi- 
tions overlaid on the best fit, with ionospheric corrections reinstated but the first epoch 
dropped. Equally weighted visibilities were used, mitigating systematic errors and ahowing 
for the first time a significant measurement of the parallax for this pulsar. 



Table 5.6. Final results (equally weighted visibilities) for PSR J1559-4438 



Parameter 



Naive fit 



Bootstrap fit 



Iiielusivc fit 



Right ascension (J2000) 
Declination {J2000) 
tia (mas yr~ ^ ) 
MS (mas yr^^) 
TT (mas) 

Nominal distance (pc) 
Nominal vt (km s~ ) 
Reduced chi— squared 
Mean epoch fit error (mas) 
Intra— epoch sys. error (mas) 
Inter— epoch sys. error (mas) 



15:59:41.526121 ± 0.000006 
-44:38:45.901849 ± 0.000020 
1.62 ± 0.08 
13.20 ± 0.02 
0.280 ± 0.048 
3570 
22S 
3.6 



15:59:35.526121 ± 0.000008 
-44:38:45.901859 ± 0.000072 
1.60 ± 0.16 
13.21 ± 0.08 
0.291 ± 0.111 
3440 
217 



15:59:41.526126 ± 0.000008 
-44:38:45.901778 ± 0.000035 
1.52 ± 0.14 
13.15 ± 0.05 
0.384 ± 0.081 
2600 
164 
1.0 
0.242 
0.259 
0.055 
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Table 5.7. Noise-added fits for PSR J1559-4438, sensitivity-weighted visibilities 



Parameter 


Da 


Db 


Dc 


Right ascension (J2000) 


15:59:41.526163 ± 0.000027 


15:59:35.526086 ± 0.000025 


15: 59:41.526025 ± 0.000041 


Declination (J2000) 


-44:38:45.902054 ± 0.000154 


-44:38:45.902131 ± 0.000163 


-44:38:45.901760 ± 0.000258 


fict (mas yr ~ "'" ) 


2.03 ± 0.46 


1.44 ± 0.59 


2.83 ± 0.79 


fJ-S (mas yr~^) 


13.39 ± 0.24 


13.40 ± 0.21 


13.00 ± 0.44 


TT (mas) 


0.109 ± 0.278 


0.270 ± 0.373 


0.621 ± 0.579 


Nominal distance (pc) 


9200 


3700 


1610 


Nominal vt (km s~ ) 


590 


237 


102 


Reduced chi— squared 


1.0 


1.0 


O.S 


Mean epoch fit error(mas) 


0.441 


0.904 


1.308 


Intra— epoch sys. error (mas) 


0.504 


1.412 


0.863 


Inter— epoch sys. error (mas) 


0.769 


0.315 


0.0 


Average single— epoch SNR 


33 


17 


9 



Table 5.8. Noise-added fits for PSR J1559-4438, equally weij 


ghted visibilities 


Parameter 


Da 


Db 


Dc 


Right ascension (J2000) 


15:59:41.526142 ± 0.000026 


15:59:35.526115 ± 0.000042 


15:59:41.526268 ± 0.000084 


Declination (J2000) 


-44:38:45.901745 ± 0.000076 


-44:38:45.902220 ± 0.000234 


-44:38:45.901798 ± 0.000420 


lia (mas yr~ "'" ) 


0.91 ± 0.38 


0.22 ± 0.99 


-0.04 ± 1.20 


fig (mas yr^^) 


13.18 ± 0.10 


13.56 ± 0.32 


12.65 ± 0.65 


TT (mas) 


0.544 ± 0.235 


0.760 ± 0.597 


2.092 ± 0.787 


Nominal distance (pc) 


1840 


1320 


480 


Nominal vt (km s""'") 


115 


85 


29 


Reduced chi— squared 


0.8 


1.0 


0.9 


Mean epoch fit error(mas) 


0.618 


0.937 


1.853 


Intra-epoch sys. error (mas) 


0.791 


1.433 


1.631 


Inter— epoch sys. error (mas) 


0.000 


0.380 


0.000 


Average single— epoch SNR, 


20 


10 


5 



Tables 15.71 and 15.81 show that while the equally weighted dataset performs better for 
Da, when the average epoch SNR is still high, its performance rapidly deteriorates as the 
epoch SNR decreases. In Dc, the pulsar was not detected in several epochs using equally- 
weighted data. The reduction in performance is less marked for the weighted datasets, 
although they are clearly still affected by systematic errors. However, if the pulsar was 
closer and the parallax larger, these systematics would be less dominant, and weighted 
datasets would allow measurement of a parallax when equally weighted datasets may be 
overwhelmed by thermal errors, to the point of not detecting the pulsar in a single epoch. 

As noted in Section 15.5.31 the use of weighted visibility points would always be opti- 
mal if the weights could include an estimate of the systematic error contribution to that 
visibility. In the absence of such an estimate, the results obtained in this investigation 
lead to the guideline that for the LBA (with typical observing conditions and calibrator 
throws), the transition region from systematic to thermal error dominated astrometry oc- 
curs when the single-epoch detection SNR falls to approximately 10 for equally-weighted 
visibilities. This is shown by the similarity of result quality for Db, where the average 
epoch SNR was approximately 10 for the equally-weighted visibilities. Alternatively, both 
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weighting regimes can be used and average total single-epoch error (formal + systematic) 
can be compared to estimate the optimal weighting scheme. Again, this is borne out in 
the simulated datasets, where a transition from systematic errors dominating to epoch fit 
errors dominating is seen as more noise is added. These guidelines have been applied to 
the results obtained in Chapter [6l 



5.7 Contributions to systematic error 



The major contributions to systematic error in VLBI astrometry include geometric model 
errors (station/source position, Earth Orientation Parameters), residual ionospheric and 
tropospheric errors, variable phase reference source structure, and image wander due to re- 
fractive scintillation. Of these, at 1650 MHz residual ionospheric errors wou ld be expected 



to do minate, despite the a priori ionospheric calibration employed (e.g. iBrisken et al 



2002). 



It is difficult to estimate the magnitude of the potential correlation between epochs of 
the residual ionospheric errors, due to the model-dependence of the ionospheric correction. 
Whilst some correlation is likely, it is unlikely to dominate, due to the significant vari- 
ability inherent in the ionosphere. Residual tropospheric errors also have the potential for 
some correlation from epoch to epoch, for example due to imperfect modeling of seasonal 
atmospheric variations. Residual tropospheric errors, however, should be considerably 
smaller than the ionospheric errors. 

Geometric model errors cause relative astrometric errors which increase with calibrator 
throw. Earth Orientation Parameters (EOPs) are well d etermined by geode tic observations 



and make minimal contributions to astrometric errors (Pradeletal 



200a). Similarly, the 



mean position of well-determined calibrators makes a minimal contribution. Some LBA 
stations, however, have position uncertainties of several cm, which could make a several 
hundred /uas contribution to systematic error. The magnitude of the error depends on 
the source declination and the calibrator-target separation. Given similar uv coverage at 
different epochs, however, the offset will be largely constant with time and is absorbed into 
the reference position of the target. Future planned geodetic observations will continue 
to improve LBA station positions, and reduce this systematic contribution. Additionally, 
as noted in Section 15.3.21 small station position errors can be corrected post-correlation, 
which offers the potential to further improve previous position fits. 

Refractive image wander is caused by large-scale fluctuations in the ISM, and can be 
estimat e d bas ed on the strength of the pulsar scattering and the scattering disk size (e.g. 



Rickett 



1990). For strong scattering, which includes PSR J1559-4438, the image wander 
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is less tha n the scatterin g disk size, if Kolmogorov turbulence is assumed for the scattering 



material (Rickett 



199d ). Thus, since Table [531 shows that scattering disk of PSR J1559- 



4438 is estimated to be only 133 /xas at 1650 MHz, the maximum refractive image wander 
is <^ 100 ;uas, and can be discounted as a source of systematic error. Table [5^ shows that 
refractive scintillation is unlikely to be significant for any of the pulsars targeted in this 
thesis. 

The variability of calibrator structure with time depends on the source chosen, but all 
compact extragalactic ra dio source s are e xpected to show some variability, with typical 
RMS values of 100 /ias (jFomalontl . l2005l ). This image wander may be correlated from 
epoch to epoch over short time periods and absorbed into proper motion fits, but over 
long times (which could be longer than as astrometric observing program), the mean 
apparent position will be constant. Some sources, such as B2201+315 and B1739+522 
(an "other" and "candidate" member of the ICRF, respectively), have shown ap parent 
motion of many hundreds of mas per year, which can persist for up to several years (jTitovl . 



20071 ). As detailed long-term information is not available for any of the reference sources 



used in this work, the impact on proper motion due to source variability is difficult to 
quantify, but the stability of the fitted models used for amplitude refinement (as discussed 
in Section I5.3.3P suggests that any reference source variability is likely to be small. As 
shown in the following chapter, formal errors for proper motion measurements are only 
accurate to several hundred /ias per year for most pulsars, and so the contribution of 
reference source variability is likely to be minimal in most cases. 

The magnitude of observed systematic errors, and the correlation of systematic errors 
with calibrator structure and angular separation from the target, are investigated further 
with the full pulsar set in Section [6.41 and used to estimate the astrometric accuracy which 
can be obtained with present and future VLBI arrays. 



6 

ASTROMETRIC RESULTS AND 

INTERPRETATION 



The data obtained for each pulsar were processed using the pipeline described in Chapter [5l 
applying the guidelines developed to minimise the total astrometric error budget. Table [6?!] 
shows the fitted parameters for each target pulsar, as well as estimates of mean thermal 
and systematic errors. Unless otherwise stated, all errors are la. 

The implications of the measured values are discussed for each pulsar in turn in Sec- 
tions [6T] aiidj6]2]_befow_j_with the analysis for PSR J0437-4715 closely following that pre- 
sented in 



De 



ler et al 



(1200811. the analysis for PSR J0737-3039A/B closely following that 



presented in lDeller et al.l (l2009al). and the analysis for PSR J1559-4438 closely following 
that presented in iDeller et al.l (j2009bl ) . In Sections 16.31 and 16.41 the revised distance and 
velocity measurements for all pulsars are used to estimate the accuracy of pulsar distance 
and velocity models, and predict the astrometric accuracy that will be achievable with 
future VLBI observations. 



6.1 Binary millisecond pulsars 



As discussed in Section 12.2.51 a variety of fascinating physics can be probed using binary 
millisecond pulsars, such as tests of GR and neutron star formation events. The three 
binary pulsars targeted in this program were chosen because uncertainty regarding their 
true distance was limiting the science which could otherwise be achieved with the system. 
While significant parallaxes were observed for PSR J0437-4715 (120o-) and PSR J0737- 
3039A/B (6(t), no parallax could be detected for PSR J2145-0750, which had previously 
been unsuccessfully targeted with the VLB A (iBrisken et al.l . l2002l ). The implications of 
the distance and velocity measurements for each pulsar are discussed below. 
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Table 6.1. Astrometric fits for all target pulsars 



Parameter 






PSR JOlOS-1431 


PSR J0437-4715 


PSR J0630-2834 


PSR J0737-3039 


Right ascension {J20aO) 






01:08:08.347016 ± 0.000088 


04:37:15.883250 ± 0.000003 


06:30:49.404393 ± 0.000043 


07:37:51.248419 ± 0.000026 


Declination (,12000) 






-14:31:50.187139 ± 0.001069 


-47:15:09.031863 ± 0.000037 


-28:34:42.778813 ± 0.000372 


-30:39:40.714310 ± 0.000099 


fj,a (mas yr ) 






75.05 ± 2.26 


121.679 ± 0.05 


-46.30 ± 0.99 


-3.82 ± 0.62 


M5 (mas yr~^) 






-152.54 ± 1.65 


-71.820 ± 0.09 


21.26 ± 0.52 


2.13 ± 0.23 


TT (mas) 






4.170 ± 1.421 


6.396 ± 0.054 


3.009 ± 0.409 


0.872 ± 0.143 


Distance (pc) 
Vt (km s"-"-) 
Visibility weighting 






Sensitivity 


156. 3+J^ 

I04.7t;;° 

Equal 


332^11 
Equal 


24+" 
Sensitivity 


Image weighting 






Natural 


Uniform 


Uniform 


Natural 


Average epoch mean fit error 


mas) 




1.232 


0.059 


0.765 


0.747 


Average intra— epoch systematic error 


(mas) 


2.477 


0.068 


0.839 


0.939 


Average inter— epoch systemat 


c error 


'mas) 


4.310 


0.103 


1.205 


0.0 


Average single— epoch SNR 










15 




Parameter 






PSR ,11559-4438 


PSR .12048-1616 


PSR J2144-3933 


PSR J2145-0750* 


Right ascension {J2000) 






15:59:41.526126 ± 0.000008 


20:48:35.640637 ± 0.000040 


12:44:12.060404 ± 0.000045 


21:45:50.461901 ± 0.000098 


Declination (J2000) 






-44:38:45.901778 ± 0.000035 


-16:16:44.553501 ± 0.000147 


-39:33:56.885041 ± 0.000316 


-07:50:18.462388 ± 0.000558 


flex (mas yr~ ""^ ) 






1.52 ± 0.14 


114.24 ± 0.52 


-57.89 ± 0.88 


-15.43 ± 2.07 


fiS (mas yr~^) 






13.15 ± 0.05 


-4.03 ± 0.24 


-155.90 ± 0.54 


-7.67 ± 0.81 


TT (mas) 






0.384 ± O.OSl 


1.712 ± 0.909 


6.051 ± 0.560 




Distance (pc) 
Vf (km s~ ) 
Visibility weighting 






2600+5|0 
Equal 


,,7-1-362 
Equal 


165±1I 

130t\t 
Equal 


Sensitivity 


Image weighting 






Uniform 


Uniform 


Natural 


Natural 


Average epoch mean fit error 


mas) 




0.242 


0.517 


1.025 


2.136 


Average intra— epoch systematic error 


(mas) 


0.259 


1.282 


1.450 


0.0 


Average inter— epoch systemat 


c error 


'mas) 


0.055 


0.106 


0.875 


0.0 


Average single— epoch SNR 






50 


23 


10 


8 



Based on two detections, with parallax fixed at 0. Variation of the parallax value between and 2 mas results in less than 100 fias yr ^ difference in derived proper motion. 
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Fitted pulsar patti - 
Fitted epocli point 
Position moasuremenl - 
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en 
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Time (l«JD) Time (MJD) 

Figure 6.1 Motion of PSR J0437-4715, with measured positions overlaid on the best fit. 
Clockwise from top left: Motion in declination vs right ascension; as before but with 
proper motion subtracted; right ascension vs time with proper motion subtracted; and 
declination vs time with proper motion subtracted. 



6.1.1 PSR J0437-4715 



PSR J0437-4715 is the brig htest and nearest observed mil 



isecond pulsar, and ha s also 



et (IKargaltsev 



Johnston et al 



et al 



200J), and 



1993)- The high 



been st udied in the optic al (jBell et al.l . Il993l ). ultravio! 
X-ray (jZavlin et al.l . 12002 ) bands since its discovery by 
rotational stability and close proximity of this pulsar-white dwarf binary system make it 
an excellent probe of General Relativity ( GR) and alternate forms o f gravitational theories. 
The measurement of its Shapiro delay by van Straten et al] ( 2001 ) is one such test which 
has shown consistency with GR predictions. The search for the low fr equency stochastic 
gravitational wave background (GWB) using pulsar timing arrays (e.g. iJenet et al.l . l2005l ) 



is another test of GR which is facilitated in part by timing of PSR J0437-4715. 



VLBI results 



Figure ETT] shows the fitted and measured positions for PSR J0437~4715, and the full set 
of fitted parameters is shown in Table [UTTl The distance of 156.3 it 1.3pc derived for PSR 
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J0437-4715 from these observations is the most accurate distance measurement (in both 
absolute and fractional distance) for a pulsar to date and approaches the most accurate 
distance meas u remen ts made of objects outside the Solar System (T Tauri, 147.6 it 0.6 pc: 



Loinard et al 



20071 ). Previo usly, the high e st-pr ecision VLBI pulsar distance determi- 

(j2002l ). who measured the distance of PSR 



Brisken et al 



nations were those made by 
J0953+0755 to an accuracy of 5 pc, along with 8 other Northern Hemisphere pulsars. The 
two previous paralla x measurements made using a southern array, of PSR J0835-4510 



(Dodson et al 



20031 ^ and PSR J 1456-6843 (jBailes et al 



1990j), had lo" distance errors of 



19 and 70 pc respectively. 

The distance to PSR J0437-4715 estimated from the TC93 model (140 pc) and the 
NE2001 model (189 pc) is within ~20% of the measu red distance of 156 . 3 ifc 1.3 pc in both 



cases . Previous timing par allaxes of PSR J0437-4715 (jHotan et al 



2006 



van Straten et al 



2001 



Sandhu et al 



19971 ) have yielded considerable variation in measured p arallax, the 



cause of which is ascribed to the inaccuracy of earlier Solar System ephemerides (jVerbiest et al, 



20081 ). These observations, which are independent of the Solar System ephemeris, provide 



confirmation that the earlier measurements of timing parallax were inaccurate. 

Previous measurements of the transverse speed of PSR J0437-4715 using scintillation 

(1199811 measured a 



Johnston et al 



observations have also provided widely varying results. 
scintil lation speed of 170 km s~^ assuming a distance of 180 pc, while lGothoskar and Gupta 
( 200ol ) measured 231km s^^, also assuming a distance of 180 pc. Neither of these values 



are consistent with the well-determined VLBI value of 104.7ib 1.0 km s~^. As discussed in 
Section 16.3. H errors in scintillation speed estimates of this magnitude are not uncommon, 
and may often be due to the simplifying assumptions made about the scattering material. 



Comparison to timing astrometry 



To compare the VLBI and timing positions, the timing data of lVerbiest et al.l (120081 ) have 
been re-fitted to obtain the position of PSR J0437-4715 at the VLBI reference epoch (MJD 
54100). Table [62] shows that the timing and VLBI positions at MJD 54100 differ by over 
two mas, many times the formal errors shown. However, the formal errors are negligible 
compared to th e uncertainty in the VLBI phase reference calibrator position (0.8 mas; 



Ma et al.l . Il998l ) , and the potential constant offsets due to phase-re ferencing errors suc h 



as station position errors (known to exist at the cm level for the LBA: lDeller et al.l . l2009bl ). 
Discrepancies between interferometri c and timing positio ns of even larger magnitudes have 
been found using the DE200 frame (jBartel et al.l . Il996l ). and differences at the mas level 
still exist for the position PSR J0437-4715 calculated using the newer DE414 solar system 
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Table 6.2. Fitted VLBI results for PSR J0437-4715 and comparative timing values 
(positions re-referenced to the VLBI proper motion epoch) 



Parameter 



Fitted value and error 



VerbiestetaL f2008) timing values 



Right Ascension (J2000) 
Declination {J2000) 
(ia (mas/yr) 
(IS (mas/yr) 
Parallax 7r (mas) 
Distance (pc) 

Transverse velocity (km/s) 

Reduced chi— squared 

Average epoch mean fit error (mas) 

Average intra— epoch systematic error (mas) 

Average inter— epoch systematic error (mas) 

Reference epoch for proper motion (MJD) 



04''37'"15. 883250" ± 0.000003 
-47°15'09. 031863" ± 0.000037 
121.679 ± 0.052 
-71.820 ± 0.086 
6.396 ± 0.054 
156.3 ± 1.3 
104.71 ± 0.95 
1.0 
0.059 
0.068 
0.103 
54100.0 



04'>37"> 15.883185= ± 0.000006 
-47° 15'09. 034033" ± 0.000070 
121.453 ± 0.010 
-71.457 ± 0.012 
6.65 ± 0.51 
157.0 ± 2.4* 
104.9 ± 1.6* 



Derived from the kinematic distance obtained from , not the less precise parallax values 



ephemeris, as compared to the DE405 ephemeris used bv lVerbiest et al.l (j2008l ). Thus, it is 
concluded that the VLBI and timing position difference is consistent with the uncertainty 
in the calibrator position and the offset between the solar system frame and the ICRF. 

The parallax value obtained from VLBI is consistent with that derived from timing, 
and yields a distance which is consistent with the kinematic distance of 157.0 it 2.4 pc 
derived from the orbital period derivative -Pb- However, the VLBI parallax measurement 
is an order of magnitude more precise than the timing measurement, and yields a distance 
which is a factor of two more precise than the kinematic distance. 

Finally, the values obtained for proper motion from VLBI differ by ~ 4cr in both right 
ascension and declination from the values estimated from pulsar timing. A likely cause 
for this discrepancy is small changes in the centroid position of the phase reference source 
due to intrinsic source variability, which would be absorbed into the astrometric fit. If 
the centroid position change is roughly linear over the timescale pr obed, the effects would 
largely be absorbed into the measured proper motion. iTitovl (j2007l ) show that some ICRF 
sources exhibit apparent proper motions of hundreds of //as yr^^ for periods of several 
years. Whilst the calibrator used here (J0439-4522) has not p r eviou sly been the target of 
detailed variability studies such as those undertaken by iTitovl ((20071), it is possible to use 
archival observations to show that centroid position changes are a plausible explanation 
for the proper motion discrepancy, as shown below. 

First, the VLBI data obtained during this thesis on J0439-4522 can be inspected di- 
rectly for evidence of variability. The model used for J0439-4522 in this thesis consisted 
of a ~ 1 mas FWHM Gaussian, with two additional delta components within 5 mas of 
strength 2% and 0.2% of peak flux. There is no gross evidence of variability over the four 
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observational epochs, as the width and positions of the primary and secondary compo- 
nents remained constant to within ~ 200 /ias, and secondary fluxes remained constant to 
~lmJy. However, as the source is only barely resolved (the beamsize is ~3mas) posi- 
tional variability at the ~ 100 //as level would be difficult to detect. An indirect indicator 
of the past variability of the source can instead be obtained from calibrator measurements 
usme; the ATCifl. These show that the total flux density on ATCA baselines has declined 
by a factor of three in a near-linear manner over a four year period from 2004 to 2007. 
On the shortest timescales probed by this archival data (several months) departures from 
the linearly decreasing trend on the order of 10% of the total flux are seen. 

The commensurate ATCA observations taken simultaneously with VLBI observations 
during this thesis are unsuited to accurate flux monitoring, due to the compact array 
configurations used and the lack of ATCA flux calibrator scans, which are unnecessary 
for VLBI. Thus, the archival ATCA observations provide the best available means for 
estimating source variability. The steady decrease in flux level is consistent with a near- 
linear change in the source centroid, which would be absorbed into the proper motion 
estimate. Smaller non-linear centroid changes are also possible, which would affect not 
only the proper motion, but also the parallax estimated from the VLBI measurements. 
However, the limited amount of short-term variability seen is consistent with a mostly 
linear change in position centroid. Thus, the archival ATCA data supports the theory 
that centroid changes in the phase reference calibrator are responsible for the difference 
between the VLBI and timing measurements of proper motion, and indicates that the 
effect of the reference source variability on parallax is likely to be considerably smaller. 

The higher proper motion precision obtained with the timing data (a factor of 5-7 
times better than the VLBI results) reflects the fact that the timing data spans a time 
baseline 7 times longer than the VLBI dataset. 



Limits on anomalous accelerations 

The newly measured parallax of vr = 6.396 it 0.054 mas allows an improved measurement 
of any anomalous acceleration of either the Solar System or PSR J0437-4715. Specifically, 
the apparent acceleration due to time variability of Newton's gravitational constant G and 
the mass of an undetected trans-Neptunian planet near the line of sight to the pulsar can 
be limited. 

As first described by iDamour and Tayloii (jl991 



a precise measurement of a binary 
pulsar's orbital period derivative, Pb) can be used to constrain a variation of the gravita- 



^http://www.narrabri.atnf.csiro.au/cgi-bin/Calibrators/calfhis.cgi?source=0437-454&band=3cm 
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tional constant, G. However, as lBell and Bailed (j 19961 ) pointed out, for PSR J0437— 4715 
a precise di s tance needs to be known in order to correct Pb for the Shklovskii acceleration 



(Shklovskii. 



1970 ) caused by it s proper motion. Thi s analysis has been performed, based 

([200a), whose limit was dominated by the 



Verbiest et al, 



exclusively on timing data, by 
uncertainty in their parallax measurement. This VLBI parallax value improves their limit 
by a factor of nearly 10 down to G/G = (-5 ± 26) x 10"^^ yr"^ at 95% certainty. This 
value compares well to the most stringent lim it currently pub l ished : (4 it 9) x 10 



-13 



Williams et al. 



200 



3). 



yr 



Since vr and Pb are 



derived through Lunar Laser Ranging (LLR; 
now both measured to similar precision, both measurements will have to be improved for a 
further significant increase in G sensitivity. Additional VLBI observations and continued 
timing could see this limit improve upon the existing LLR limit early in the next decade. 

An alternative source of anomalous acceleration is hea vy planets in a wide orbit aroun d 
the Sun or the pulsar. Building upon the initial analysis of lZakamska and Tremaind (120051) 



this p arallax measurement can be combined with the timing results from 



Verbiest et al 



(|2008l ) to derive the following result: ap,/c = (3 ± 16) x 10"^" s"^ at the 2a level. This 
improves the limit published in I Verbiest et al. (2008) by an order of magnitude and makes 
PSR J0437-4715 a more sen sitive Solar System acceleromet er than PSR J1713+0747, the 
most precise pulsar listed bv lZakamska and Tremaind (|2005l ). From this, the limit for 50% 
of the skj^ can be calculated as: |a0,5o%sky/c| < 3.9 x 10~^^s~"'^ (95% certainty). This 
acceleration limit can be used to exclude massive bodies within a given radius of the Sun; 
for example, at Kuiper-belt radii (50 AU) it excludes a planet more massive than Uranus 
over 50% of the sky, while Jupiter-mass planets are excluded within 226 AU over 50% of 
the sky. 



Impact of the stochastic GWB 

Using tools recently developed by 



Hobbsetal 



([2003), the effect of a stochastic GWB 
on the observed value of Pb from pulsar timing has been simulated. The c haracteristic 



strain spectrum of the GWB was set to that of the best published GWB limit (jJenet et al 
20061 ). with dimensionless amplitude A = 1.1 x 10^^^ at a period of one year and power-law 



dependence on frequency with exponent a = —2/3 (as predicted for a GWB dominated 
by black hole-black hole mergers). The simulated GWB causes the kinematic distance 
to be inconsistent with the VLBI parallax distance at the 2a level in ~ 50% of trials. 
Thus, although these observations cannot improve upon the present GWB limit, they 
are consistent with it. Simulations with a GWB with amplitude of 1.1 x 10"^'^ show 



within 60° of the hne of sight towards and away from PSR J0437-4715 
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inconsistencies between the kinematic and VLBI distances at the 2a level in 95% of trials, 
providing an independent exclusion of a GWB with an amplitude at or above this value. 

It is also interesting to note that the precise limit on G presented above would be 
impossible in a Universe with a strong GWB. In the simulations with GWB amplitude 
1.1 X 10^^'^, the observed G value is inconsistent with in 99% of cases, merely due to the 
GWB-induced corruption of the timing measurements. Thus, the stochastic GWB must 
eventually limit the accuracy of measurements of G in the fashion outlined here. 



6.1.2 PSR J0737-3039A/B 



The double pulsar system PSR J0737-3039A/B (jBurgav et al 



2003 



Lvne et al.1 . 1200J) 



is one of only eight known double neutron star (DNS) systems, and the only system in 
which both neutron st ars are visible a.s pul sars. The discovery of the mildly recycled "A" 
pulsar was reported by I Bur gay et al.l (120031 ). but emission from the companion "B" pulsar 



was not ob served at the time , and it was not discovered until follow-up observations some 



time later ( Lvne et al 



20041 ). PSR J0737-3039A/B is the most relativistic known binary 



system, with an orbital period of only 2.5 hours and a coalescence time (due to orbital 
energy loss to gravitational radiation) of 85 Myr. PSR J0737-3039 A/B provides the best 
current tests of General Relativity (GR), with lKramer et al.l (j2006l ) recently showing that 



measurements of the post-Keplerian parameter "s" agreed with the GR prediction to 
within 0.05%. The discovery of PSR J0737-3039A also led to a marked upwa rd revision 
in the estimated Galactic merger rate of DNS systems ( Kalogera et al.l . 2004), although 
merger rate estimates are still very uncertain due to the poor constraints available on the 
characteristics of the DNS population. 

One o f the major ou t stand ing uncertainties surrounding PSR J0737-3039A/B is its 



distance. iKramer et al 



(|2006l ) detect a marginally significant timing parallax of 3 it 2 
mas, but this is no more accurate than what is typically assumed for distances estimated 
from dispersion measures in the NE2001 distance model. The NE2001 DM distance 
estimate for PSR J0737-3039A/B is 480 pc, while th e earlier TC93 mode l estimates the 



Kramer et al 



200e 



Burgav et al 



dista nce to be 570 pc. Many previous studies (e.g. 
20031 ) have assumed the rounded distance of 500 pc when calculating luminosity and 
kinematic values. As discussed below, the accurate distance and velocity information for 
this system presented here enables tests of GR to proceed to greater precision in the future, 
as well as providing insight into its formation and x-ray emission. 
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VLBI results 

As shown earlier in Table I5.H PSR J0737-3039A/B was clearly the most challengin, 
target undertaken in this program - with an equivalent gated 1600 MHz flux densitjl 
of 4 mJy, it is almost as faint as PSR J0108-1431, but with a DM distance four times 
greater, a factor of four greater precision was required to attain an equivalent parallax 
accuracy. Whilst a very nearby (angular separation ~ 20 arcminutes) calibrator was 
available (B0736-303), Figure 15.21 shows that it presented a particular challenge, as it 
was weak (correlated flux density ~30mJy on 1000 km baselines), possessed complicated 
structure, and its position was not known precisely in the ICRF. As noted in Section [5.41 
the position of B0736~303 was refined using a phase-referenced VLBI observation, and 
finally determined to approximately five mas accuracy by coi nparison of the VLBI position 



for PSR J0737-3039A/B with its published timing position (jKramer et al.l . l2006l ). Whilst 
this is a relatively large positional uncertainty, the short calibrator throw means that the 
resultant differential errors are small. 

The complex structure and positional uncertainty in the phase reference source in- 
troduces an extra element of model-dependence to the results obtained for PSR J0737~ 
3039 A/B. The final result obtained for parallax differed by 2a from the initial result, 
which was obtained when the calibrator position was in error by approximately 50 mas. 
After the calibrator position was updated, several different calibrator models were pro- 
duced (model A, which used Gaussian model components; model B, which used standard 
CLEAN components; and model C, which used a combination of the two) and the results 
compared. In each case, the model was held constant over all epochs. The clean images 
of the combined dataset for B0736-303 (visibilities from all epochs concatenated into a 
single file) are shown in Figure [6?2] for models A, B and C (which is the model shown in 
Figure [521 but is shown again here for clarity). 

As shown in Figure [6^21 model C gives a significantly better fit to the data, resulting 
in lower image residuals than model A or model B. The data reduction for PSR J0737- 
3039A/B was repeated three times using models A, B and C, and the resultant position 
datasets were fit for parallax and proper motion. The results of this analysis showed 
that model C gave a significantly better fit (vr = 0.87 it 0.14 mas) than model A (vr = 
1.06 ± 0.18 mas) or B (tt = 0.76 ± 0.20 mas). 

Model C was thus clearly the best source model, and was chosen for the final results 
presented in Table 16. 1[ Since any residual model-dependent errors are likely to be small 
(and difficult to quantify in any case, but on the order of the formal errors or smaller) 



Since pulsar A dominates the radio emission, gating was performed using the pulsar A ephemeris only 
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Figure 6.2 Images of B0736-303, using Gaussian model components only (model A, top 
left), CLEAN components only (model B, top right) and a combination of Gaussian and 
CLEAN components (model C, bottom left). Images contours are at 1, 2, 4, 8, 16, 32 and 
64% of peak flux density in each case. Model C clearly gives the best fit to the data. 
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the error given for PSR J0737-3039A/B in Table EH] and the discussion below is simply 
the "inclusive" astrometric fit error described in Section 15.41 and does not include the 
uncertainty due to calibrator structure discussed here. However, it is worth noting that 
the some uncertainty due to the calibrator source structure remains, and if the model 
could be improved by additional observations with better uv coverage in the future then 
the astrometric data could be re-reduced, with the potential for a slightly improved fit. 

Figure shows the final fitted and measured positions for PSR J0737-3039A/B. The 
measured parallax for PSR J0737-3039A/B (0.872 it 0.143 mas) corresponds to a distance 
of 1150l|gQ pc. The successful measurement of parallax and proper motion by this program 
represents the first significant parallax and proper motion determination for a DNS system. 
At 95% confidence, the lower limit for the distance is 860 pc - still considerably in excess 
of the predicted DAd distance for the system. The revised value of electron density along 
the line of sight to PSR J0737-3039A/B is 0.043 cm~^, which suggests the influence of 
the Gum nebula along this sightline is less than originally thought. 

The measured transverse velocity for PSR J0737-3039A/B is 24^g km s~^. However, 
whe n the best estimates of peculiar motion of the Solar System and Galactic rotation 
(e.g. iMignardl . 12000) are subtracted in order to obtain a velocity in the local standard 
of rest, the resultant transverse velocity is 91*^3 km s~^. This is well within the range 



of transverse velocities expecte d for the massive stars w 



rich are progenitors for NS-NS 



w 


rich c 




1965) 



binaries (~ 20km s~ ; see e.g. iFeast and Shuttleworthl . Il965l ). and such a low value is 
extremely unlikely unless both kicks imparted on the system by the supernova events 
which formed the pulsars were extremely small. If the PSR J0737-3039A/B system does 
have a large velocity, it must be predominantly in the radial direction - a possibility 
discussed in detail below where possible formation models for the system are considered. 



Comparison to timing astrometry 



The parallax measurement made throu gh these VLBI obser vations is consistent with the 
timing value of 3 it 2 mas obtained by Kramer et al.l ( 2006 ). but it is considerably more 
precise. The values obtained for proper motion (^^ = —3.82 it 0.62 mas yr~^, fis = 2.13 it 
0.23 mas yr~^) are also consistent with those obtained from timing (//„ = —3.3 it 0.4 mas 

yr^^, ^5 = 2.6 ± 0.5 mas yr~^). 

Th e measured transverse velocity of 24~^q km s"'^ is over twice that obtained by[ 



Kramer et al 



(|2006l ). almost entirely due to the distance revision. However, this value is sti ll considerably 
lower than those obtained through measurements of scintillation velocity bvlRansorn et al 
( 20041 ) . who obtained 141 it 8km s~^, and the revised value of 



Coles et al 



3)^ 



20051) ■ who 
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Figure 6.3 Motion of PSR J0737-3039A/B, with measured positions overlaid on the best 
fit. From top: Motion in decHnation vs right ascension; right ascension vs time with 
proper motion subtracted; and declination vs time with proper motion subtracted. 
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obtained 66ibl5km s ^ . Both of these studies assumed a distance to PSR J07 37-3039A/B 



close to the NE2001 value (570 pc). Thus, as noted by iKramer et alJ (j2006l ). it is likely 
that the usual scintillation velocity assumption of a thin scattering screen with an isotropic 
turbulence distribution is not appropriate in this case. Accounting for anisotropy in the 
ISM turbulence gives a measurement of scintillation velocity that is lower still, and more 
consistent with the VLBI transverse velocity (W. Coles et al., in preparation). 

Since the timing position of PSR J0737-3039A/B was used to refine the VLBI cali- 
brator position, no meaningful comparison can be made between the timing and VLBI 
positions. When corrected for proper motion to the same epoch, the timing and VLBI 
positions have been aligned to approximately one mas, which is the the level of errors due 
to proper motion uncertainty. 

Finally, it is appropriate to note that the continuation of timing observations of PSR 
J0737-3039A/B should see a determination of parallax and proper motion to equal or 
greater precision than these VLBI results within the next decade. 



Implications for tests of GR 

One of the tests of GR possible with the PSR J0737-3039A/B system is the rate of 
change of orbital period (Pb) due to the loss of energy to gravitational radiation. In 
order to perform this test, the observed Ph must b e measured as accurately as possible 



Kramer et al 



200^ ). and contributing factors to Pb other 



(the current significance is 70a ; 
than GR must be estimated as accurately as possible. For PSR J0737~3039A/B, the major 
contributing factors are differential Galactic rotation (-P™'), acceleration towards the plane 
of the Galaxy (P^), and the appar ent acceleration caused by the transverse motion of the 



Shklovskii 



1970 ). These will collectively be referred to 



system (the Shklovskii effect ib'^'^' 
as galactic and kinematic contributions (-Pb*^)- At the newly calculated distance of 1150 
pc, the magnitude of thes e effects can be calculated using Equations 2.12 and 2.28 from 
Damour and Taylor! (Il99ll ) as: 



prot 

Ph 
pShk 



^ X ( cos Z + 



sin 6 



c 



sin^ I + /32 



(6.1) 
(6.2) 
(6.3) 



where I, b and z are the pulsar's Galactic longitude, latitude and height respectively, d and 
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/i are pulsar distance and proper motion, and Rq are the Galacti c speed and radius of th e 
Solar System (taken to be 7.5 kpc and 195 km s^^ respectively; e.g. Dias and LeoinJ^OOS ) 



K^, is the vertical gravitational potential of the Galaxy (taken from iHolmberg and Flvnn 



2OO4I as 0.45 km^ s"^ pc"^ at the height of PSR J0737-3039A/B) and c is the speed of 
light. The dominant uncertainties are in d (16%) and ^ (15%), and i?o, vq and Kz, whose 
errors are typically estimated at ~10%. 

Equations O - E3 give Pl°^/Ph = (-4.3 ± 0.7) x lO^^Og-i^ pzjp^ ^ (3.8 ± 0.8) x 
10~2i s-i, and P^'^^/Ph = (5.3 ± 1.8) x lO'^o g-^. Combining these terms gives P^^ / P^ = 
(1.5 ifc 2.1) X 10^^° s "^, and multiplying by the observed orbital period Pb = 8834.5 s 



(Kramer et al. 



derivative: P^ 



2006 ) yields the net effect of these terms on the observed orbital period 
= (1.3 ± 1.8) X 10-^6. 



These contributions to orbital period derivative are four orders of magnitude below 
the GR contribution, and two ord ers of magnitude be low the current measurement error 



obs 



(-1.252 ±0.017) X 10' 



-12. 



Kramer et al 



20061 ). Thus, with the current measure- 



ment accuracy of distance and transverse velocity, GR tests with PSR J0737-3039A/B 
using Pb can ultimately be made to the 0.01% level, althourfi of order ten y ears of contin 



Kramer et al. 



ued p recision timing will be required to reach this accurac}0. As noted by 
()2006l ). measuring Pb at this level will place stringent requirements on the class of gravi- 
tational theories which predict significant amounts of dipolar, as opposed to quadrupolar, 
gravitational radiation, exceeding the best solar system tests. Measuring the moment of 
inertia of pulsar A via a measurement of the spin-orbit coupling of the system, however, 
wo uld require a further order of magnitude improvement in the measurement precision of 
Pb (jKramer et al.l . 120061 ). In the near future, additional VLBI and/or timing measurements 
can be expected to reduce both the distance and velocity errors to PSR J0737-3039A/B 
below 10%; however, even with negligible error in these parameters, the existing accuracy 
of measurements of Rq, vq and would limit the accuracy of Pb measurements in this 
system to 0.004%. To attain the 10"^ precision necessary to measure the neutron star mo- 
ment of inertia, the constants Rq and vq must be measured to a precision approaching 1%. 
Such measurements may be possi ble within a deca de via continued optical observations 



of stars near the Galactic center (iGhez et al. 
likely to be challenging. 



20081 ). although attaining this precision is 



^since the Pb err or scales with t and the current significance of 70(t was reached in 2.5 years 
(jKramer et all . I2OO6I ) 
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Implications for formation models 



As noted in Chapter ^ most neutron stars are believed to receive large kicks at birth. 
Indeed, the discovery of the first binary pulsar PSR B1913+16 with its large orbital 
eccentricity (e = 0.617) supported this view and led to a "standard model" in which 
the newly-formed neutron stars in compact bi naries received large kicks and wer e often 



1991 



disrupted by this and the associated mass loss (iBhattacharya and van den Heuvell . 

PSR J0737-3039A/ B, however, does not support the standard model. Its low orbital 
eccentricity of just 0.08 (jKramer et al.l . 120061 ) and low transverse velocity suggest that very 
close binaries might influence the evolution of the progenitor stars and prevent large mas s 
loss and kicks d uring the supernova events (jPodsiadlowski et al.l . l2005l : IStairs et al.l . 120061 ). 



Willems et al, 



(120061 ) have argued that the observed properties of PS R J0737-3039A/B 



Piran and Shaviv 



can st ill be explained by large kicks and significant mass loss, however 
(|2005l ) suggest almost no mass loss took place beyond the binding energy. The refinement 
of the local transverse velocity of the PSR J0737-3039A/B system presented here allows 
a comment on the competing formation models for the system. 

Since, as shown above, the transverse velocity of PSR J0737-3039A/B is no greater 
than the likely peculiar velocity of its progenitor system, if the PSR J0737-3039A/B 
system received a large velocity kick, it must possess a large radial velocity. However, in a 
DNS system there are no observational methods available to determine the radial velocity. 
Due to t he accurate measurem ent of its Shapiro delay, PSR J0737-3039A/B is known to lie 
edge-on (jKramer et al.l . 120061 ) , and if the only kick provided to the system was provided by 
the loss of binding energy during the supernova explosion, the resultant 3D space velocity 
should be on the order of ^50 km s~^, estimated from the system's observed eccentricity 
and orbital velocity (e.g. iBlaauwI . Il96ll ). This velocity would be constrained to the plane 
of the orbit and simple geometry allows the estimation of the probability of observing a 
transverse velocity less than 10 km s~^ to be about one in eight, which is small, but not 
unreasonable. 

Conversely, if the double pulsar had received a large kick (e.g. 



Willems et al. 



20061 V 



then the odds of observing such a low transverse velocity become increasingly remote. 
Not only would the radial velocity have to be increasingly large, but the inclination angle 
of the system must not be altered by the kick. Hence, the transverse velocity results 



shown here reinforce those of 
of 



Kramer et al 



(|2006l ) and strongly favour the interpretation 



Piran and Shavivl (120051 ). who argue for almost no mass loss and kick in the case of 



PSR J0737-3039A/B. 

An intriguing aside to the implication of low kick velocities in PSR J0737-3039A/B- 
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like systems is the possible, albeit speculative, explanation this offers to the formation of 
the newly-discovered binary pulsar PSR J1903-I-0327, a fast, heavy, and highly recycled 
millisecond pulsar ( mass 1.8 solar ma,sses, period 2.15 ms) with a one solar mass main se- 
quence companion (IChampion et alJ . l2008l ) . The formation of this system has confounded 
existing models, since it possesses an intermediate orbital eccentricity of e = 0.44, and the 
orbit of a highly recycled pulsar such as PSR J1903-I- 0327 should have been circularised 
during the mass-transfer phase to a very high degree (lAlpar et al, 



If the progenitor system to PSR J0737-3039A/B had originally possessed a solar mass 
companion in a wide orbit as part of a triple system, then the kicks due to the PSR 
J0737-3039A/B supernovae would not have been sufficient to disrupt the triple system. 
In such a situation, PSR J0737-3039A and B would coalesce in 85 Myr due to the emission 
of gravitational radiation, creating either a heavy millisecond pulsar or a black hole that 



continues to orbit the companion star. 



van den Heuvel and Bonsema 



(|l984l ) ] proposed 



such a merger mechanism for the creation of the first millisecond pulsar shortly after its 
discovery. In this scenario, the majority of the binding energy of the relativistic binary 
will be lost in just a few seconds, and the resultant heavy millisecond pulsar or a black 
hole is left in an eccentric orbit with the remaining main sequence star. The eccentricity 
would be determined by the amount of matter lost. 

Such a system would strongly resemble PSR J 1903-1-0327, but only if the coalescing 
neutron stars formed a millisecond pulsar with a low magnetic field, and approximately one 
solar mass was lost from the coalescing system in much less than an orbital period of the 
outer star. The PSR J0737-3039A/B system possesses ample angular momentum to form 
a fast millisecond pulsar, but the enormous energy lost as the neutron stars are disrupted 
may help regenerate the magnetic field of the resultant heavy neutron star. Nonetheless, 
the low space velocity of the PSR J0737-3039A/B system implies that such a formation 
pathway is realistic, and that it could potentially explain the PSR J1903-I-0327 system. 



Implications for DNS merger rates 



The size of the Galactic population of DNS systems, and hence DNS merger rate estimates, 
depends upon a host of assumptions including the luminosity functions and beaming frac- 
tions of recycled pulsars. As already noted, the poor constraints on the characteristics of 
recycled pulsars mean t hat the merger rate e stimates remain uncertain. In many recent 
merger rate studies (e.g. iKalogera et al.l . l2004l ) the luminosity function of recycled pulsars 
(and in particular, its faint tail) has been inferred from the entire pulsar population, avoid- 
ing the small-number problem which would be inherent in attempting to model the faint 
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end of the recycled pulsar luminosity function directly. However, this approach would not 
be valid if the luminosity function of recycled pulsars differs significantly from that of the 
general pulsar population. Typically, multiple population models have been considered, 
leading to the widely varying merger rate estimates. 

The distributions of beaming fraction and pulse sha pe for recycled pulsar s are known 



to differ considerably from those for slower pulsars (e.g iKramer et al.l . Il998l ). and hence 
it is plausible that their luminosity distribution also differs markedly. The observational 
constraints on the luminosity functions for normal and recycled pulsars are presented in 



Lorimeij (|2008l ). which highlights the considerable extrapolation from the small observed 



sample of faint recycled pulsars. The revised distance presented here shows that the 
radio luminosity of PSR J0737-3039A/B is considerably greater (by a factor of five) than 
previously assumed. If this revision were to markedly influence the assumed luminosity 
distribution for recycled pulsars in DNS systems, then the assumed space density of DNS 
systems would be reduced, since there must be fewer systems in the Galaxy to match 
the number of observed systems. This would demand that estimates of the Galactic DNS 
merger rate and c orresponding de t ection rates for LIGO and Advanced LIGO such as 
those presented in iKalogera et al.l (j2004l ) also be reduced. However, it is infeasible to 
make conclusive statements on the basis of the single measurement presented here, and 
the true recycled pulsar luminosity function is likely to remain the object of intense study 
for some time to come. 



Implications for x ray production 



A final area in which this distance determination can contribute to understanding of the 
PSR J0737-3039A/B system is its x-ray emission. PSR J0737-3039A/B has been observed 



Pellizzoni et al 



intensively in x-rays ( McLaughlin et al.l . |2004| : IChatteriee et al.l . 120071 : 



Possenti et al 



200 



20081 ). with the x-ray luminosity measured at (2.3 — 2.4) x lO^'^erg s ^, 



assuming the DM distance of 500 pc. The newly measured VLBI distance of 1150 pc 
means this luminosity is revised upwards to 1.2 x 10^^ erg s~^. The maj ority of the x-ray 



emission is seen to be modulated at the spin period of the A pulsar (jChatterjee et al 



20071 ). However, there has been considerable debate over where and how the observed 
x-rays are generated, given the small binary separation and interacti on of the pulsar 
wind of the A pulsar with the magnetosphere of the B pulsar (see e.g. iLvutikovl . 12004 ) 



which provide alternate x-ray genera tion mechanis ms to the commonly assumed mag- 
netospheric and thermal origins (e.g. iBecked . l2000l ). Generation of x-rays due to such 



wind — magnetosphere interactions would be expected to show orbital modulations, but 



124 



Chapter 6. ASTROMETRIC RESULTS AND INTERPRETATION 



no su ch orbital modulation had been seen until a tentative detection by iPellizzoni et al 
([2003) of x-ray emission modulated at the period of the B pulsar over part of the orbit. 



Magnetospheric x-r ay emission is believed to occ ur in many pulsars that display non- 
thermal x-ray spectra (jBecker and Triimpen Il999l ). and there is no reason to believe 
that magnetospheric x-ray emission should not be generated in the PSR J0737-3039A/B 
system. Since pulsar A has a spin-down luminosity orders of magnitude higher than pulsar 
B (5.9 X 10^^ erg s^^, comp ared to 1.7 x 10 '^ ^ erg s~^) it would be expected to dominate any 
magnetospheric emission. Possenti et al. (200a) state that the if the x-ray luminosity of 
the PSR J0737-3039A/B system is dominated by magnetospheric er nission from pulsar A , 
the predicted x-ray luminosity (based on relation Lx cc E^'^ from iGrindlay et al.ll2002l ) 
is 5 X 10'^'^ erg s"^. Given the typical scatter of an order of magnitude about the Lx — E 
relation, they concluded that this is consistent with the then-presumed x-ray luminosity 
of 2.3 X 10^*^ erg s~^. Counting against the magnetospheric origin, on the other hand, 
was the fact that at the DM distance of 500 pc, the neutral hydrogen column density 
calculated from the absorbed fit (A'^h ^ 1-5 x lO^'^ cm~^) was higher than expected from 
the pulsar's presumed location in the Gum nebula (a nearby, large HII region which is 
modeled at a distance of ~ 500 pc with a depth of several hundred pc in the NE2001 
model). 



However, as noted by 



Possenti et al 



(120081 ). the value of A^h is consistent with the 
usual average of 10 neutral hydrogen atoms for every free electron along the line of sight. 
The revised distance estimate presented here places PSR J0737-3039A/B beyond the Gum 
nebula, making the Ah measurements consistent with expectations. The revised value of 
x-ray luminosity (1.2 x 10^^ erg s~^) is approximately a factor of two from the predicted 
value, which is no more discrepant than the original estimate. Hence, this revised distance 
measure strongly supports a power-law model of magnetospheric origin (from pulsar A) 
for the bulk of the x-ray emission from the PSR J0737-3039A/B system. 

Whilst the original spin-down luminosity to x-ray luminosity conversion efficiency of 
0.04% was consistent with expectations, the revised conversion efficiency of 0.2% implied 
by this new VLBI distance is the highest amongst observed millisecond pulsars (see e.g. 



Zavlin 



20061 ). There is considerable scatter in the observed values of x-ray conversion 



efficiency amongst millisecond pulsars, however, and so this is not taken as evidence of an 
unusual x-ray production mechanism. It is also worth noting that at the 95% confidence 
lower limit of distance (860 pc) the x-r ay conversion efficiency would be 0.12%, less than 



PSR B1937+21 and PSR J0218+4232 ( Zavlin 



20061 ^ 
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6.1.3 PSR J2145-0750 



PSR J2145-0750 is a binary p ulsar with sp i n per iod 16.05 ms and orbital p eriod 6.84 
days, which was discovered by iBailes et alj (j 19941 ) . Optical observations by 



Bell et al 



(jl995l ) tentatively ident ified the expect e d wh ite dwarf companion to PSR J2145-0750, 
which was confirmed by lLundgren et alJ (119961^. The scintillation ve l ocity of PSR J2145- 



Nicastro and Johnston 



Johnston et al 



1995), and a revised 



(11993). Scintillation 



0750 was measured to be 31 it 25 km s~^ by 
value of 51 km s ~^ (no error given) wa s publi shed by 
observations by iGothoskar and Guptal (j200d ) suggested a considerably higher transverse 
velocity of 113 km s~^, but this was based on a single observation and may have been 
biased by refractive scintillation effects. All scintillation studies used the TC93 distance 
of 500 pc for this pulsar. 



Timing astrometry using a 10 year dataset bv lLohmer et al.l (j 20041 ) measured a parallax 
of 2.0ib0.6mas, and a proper motion of 14.1ib0.4mas yr~^, corresponding to a transverse 
velocity of 33 it 9 km s~^. The implied distance of 500]'j^}5 pc is consi stent with th e DM 



based estimates of 500 pc (TC93) and 566 pc (NE2001). However, iHotan et al.l (120061 ^ 
reported no detection of parallax and an upper limit of 0.9 mas, at 95% confidence, 
albeit with a shorter (2.5 year) dataset, in addition to a proper motion measurement 
of 13.5 lb 6.0 mas yr"^. A previous att empt has been made to measure the VLBI parallax 
of this system by 



Brisken et al 



(I2OO2I ) using the VLBA, but the pulsar was not detected 
and dropped from the observing program. Accordingly, the distance to this system is still 
controversial and confirmation of a previous result is keenly sought. 

From the four observations made of PSR J2145-0750 during this observing program, 
signific ant detections were m ade on only two occasions. As was the case for the VLBA pro- 
gram of 



Brisken et al 



(|2002l ). it is believed that refractive scintillation is the major cause of 
the non-detections, as the pulsar's flux varies significantly over long timescales. Although 
measurement of the parallax of PSR J2145~0750 was thus not possible, by holding paral- 
lax fixed at zero mas it was possible to measure a proper motion of fia = —15.4 it 2.1 mas 
7.7 ifc 0.80 mas yr~^. Assuming a distance of 500 pc (both the TC93 esti- 



yr 



^5 



mate and the 



Lohmer et al 



(120041 ) timing parallax measurement), this co rresponds to a 



Lohmer et al 



transv erse velocity of 40 km s^^, 20% higher than the result given by the 
(|200J) timing. Varying the parallax between zero and two mas results in a change in 
proper motion of less than 0.1 mas yr~^, much smaller than the formal errors. As the 
proper motion and position have been derived from only two measurements, there are no 
degrees of freedom to the fit and no estimate of systematic errors, and so the result should 
be treated with some caution. The fitted and measured positions of PSR J2145~0750 are 



126 



Chapter 6. ASTROMETRIC RESULTS AND INTERPRETATION 




Figure 6.4 Motion of PSR J2145-0750, with measured positions overlaid on the best fit. 
Since only two significant measurements were made, the parallax of the pulsar was held 
fixed at zero, and the fit to position and proper motion has no free parameters. Variation 
of the parallax between zero and two mas (corresponding to a distance > 500 pc) made 
an insignificant (< 100 /ias yr~^) difference to fitted proper motion. 



shown in Figure EH 

The observed total proper mo tion of 17.2 ± 2.2 m as yr^^ is consistent with the value 
of 13.5 lb 6.0 mas yr~^ obtained by lHotan et al.l ( 2006), but only c onsist ent at the 2a level 

As these VLBI 



Lohmer et al 



(200 



with the value of 14.1 it 0.4 mas yr~^ obtained by 
results are derived from only two measurements, and there is no way of quantifying the 
systematic error, it is likely that the timing results are more reliable for this pulsar. 



6.2 Isolated pulsars 

Whilst isolated pulsars do not generally allow the tests of GR which can be made using 
binary pulsars, many are interesting for other reasons. Apart from providing more (and 
generally easier) targets for studying the accuracy of DM distance models, as discussed in 
Section 12.3.11 isolated pulsars can also be used to investigate the accuracy of scintillation 

""since PSR J2145-0750 lies close to the ecliptic, timing observations obtain large covariances between 
proper motion in right ascension and declination, and so the comparison of total proper motion is more 
useful 
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velocity measurements and the ali gnment of veloci t y and rotation vectors in the general 



pulsar population, as described by iJohnston et alJ (j2005l ). Another comparison that can 



be made is that of "kinematic age" (the time it would have taken the pulsar to reach its 
present Galactic latitude, given its present proper motion) to characteristic age, which can 
yield some insight into the birth location and/or spin period of pulsars. However, apart 
from the "technique check" pulsars PSR J1559-4438 and PSR J2048-1616, the pulsars 
included in this sample were chosen primarily for their unusual luminosity characteristics. 
In each case, an error in the assumed DM distance measurement could have potentially 
explained these unusual luminosities, and these VLBI observations offered a mechanism 
to remove this uncertainty. 

At the outset of this analysis, it is appropriate to note the importance of geometric 
considerations in determining pulsar luminosities, particularly in the radio. The extremely 
anisotropic nature of pulsar radio emission, due to the (non-uniform) beamed emission 
over a finite solid angle, makes direct comparison between pulsars problematic. For any 
given pulsar, the emission beam may sweep directly over the observer, or the edge of the 
beam may only graze the line of sight, with an obvious impact on the apparent luminosity. 
Pulsar apparent luminosities are generally quoted in mJy kpc^, where the pulsar flux is 
simply scaled by distance squared, and no attempt is made to account for beaming effects. 
Whilst beaming effects make radio luminosity analyses problematic for individual pulsars, 
studies of pulsar radio luminosity functions are nevertheless extremely important as the 
pulsar luminosity function is an essential component of pulsar simulations. 

The results obtained for each isolated pulsar observed in this program are discussed in 
turn below. 



6.2.1 PSR J0108-1431 



Tauris et al 



(jl994l ). and was postulated to be the 



PSR J0108~1431 was first reported by 
nearest observed radio pulsar - its DM of 2.38 remains the lowest of any known pulsar. 
In the TC93 Galactic electron model, PSR J0108-1431 was pr edicted to lie 130 pc from 
Earth - in the newer NE2001 model (jCordes and Lazid . l2002l ). the distance was revised 
to 184 pc, which if correct would double the apparent radio luminosity, to a value of 300 
//Jy kpc^ at 400 MHz (the lowest known value is that of PSR J0006+1834, which at a 
-DM-derived distance of 700 pc has an apparent 400 MHz luminosity of 100 ^Jy kpc^). 

Figure [631 shows the fitted and measured positions for PSR J0108-1431. As shown 
in Table WA\ the measured distance of 2^^^^ pc is consistent with the NE2001 distance, 
but inconsistent at 95% confidence with the earlier TC93 distance. It confirms that PSR 
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Figure 6.5 Motion of PSR J0108~1431, with measured positions overlaid on the best fit. 
From top: Motion in dechnation vs right ascension; right ascension vs time with proper 
motion subtracted; and dechnation vs time with proper motion subtracted. 
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J0108-1431 is more distant than originally suspected, and while its apparent radio lumi- 
nosity is still low, it is no longer remarkably so. In fact, 5 other pulsars (PSR J0006-I-1834, 
PSR J1829+2456, PSR J1918+1541, PSR J2015+2524, and PSR J2307+2225) have lower 
400 MHz luminosities than the revised value of 510 //Jy kpc^ for PSR J0108-1431. Of 
course, none of these pulsars have an independent distance measure, so it is possible that 
one or more of them may also be more distant than the model predictions, and hence 
also have greater luminosities. In addition, as noted above, beaming effects mean that 
these measured apparent luminosities may differ considerably from the true luminosities 
for these pulsars. 

This higher distance, however, does imply an even lower value of electron density than 
that predicted by the NE2001 model: 0.0099 cm~'^. This is amongst the lowest values 
amongst pulsars within one kpc of the Solar System, and is the lowest value for pulsars 
nearer than 500 pc. 

The detection of a parallax and proper motion enables for the first time a calculation 
of the Shklovskii correction for PSR J0108-1431, as shown: 



p 



c 

170^ X 240 
2.998 X 108 



(6.4) 



2—2 —1 

mas yr pc s m 



2.31 X 10~2 X 7.28 X 10"^*^ s"^ 



1.68 X 10"^^ s"^ 



0.808 X 1.68 X 10 



-17 



1.37 X 10"^^ 



(6.5) 



where D is the pulsar distance, fi is the total pulsar proper motion and c is the speed of 
light. 

Since the measured P for PSR J0108-1431 is 7.704 x lO^^'^ (jHobbsetalJ, boOJ) , the 



Shklovskii term contributes 18% of the observed rate of change of spin period for this 
pulsar. The consequent downward revision to the intrinsic spin-down luminosity of PSR 
J0108-1431 is important in the interpretation of the optical and x-ray conversion efficiency 
for PSR J0108-1431, as shown below. Additionally, the characteristic age of PSR J0108- 
1431 is revised upwards from 166 x 10^ yr to 200 x 10^ yr. 

PSR J0108-1431 was recently detected in x-rays by 



Pavlov et al 



(|2008l l. who found 



that the observed emission was compatible with a magnetospheric origin, but that a ther- 
mal component from heated polar caps could not be excluded. The spin-down luminosity 
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to x-ray luminosity conversion efficiency was 0.4% when assuming a distan ce of 130 pc 
This value is much larger than for young pulsars (typically 10~^ to 10~ 



Pavlov et al. 



20081 and references therein), but comparable to most other old, nonrecycled pulsars. At 



the best-fit distance for PSR J0108-1431, the increased x-ray luminosity coupled with the 
downward revision to spin-down luminosity increases the conversion efficiency to 1.7%. 
This would make PSR J0108-1431 the most efficient known pulsar in terms of converting 
spin-down luminosity into x-ray luminosity, as the distance to the previously most efficient 
pulsar, PSR J0630-2834, has been revised dramatically downwards by the measurements 
of this thesis. PSR J0630-2834 is discussed further in Section [6. 2. 21 below. However, at the 
95% confidence lower distance estimate (143 pc) the x-ray conversion efficiency estimate 
is little changed at 0.5%. A improvement on the parallax measurement would allow more 
useful constraints to be placed on the x-ray conversion efficiency. 



Deep optical searches for PSR J0108-1431 have been made using the VLT (iMignani et al 



20031 ) and while the puls ar was not detected in the original observation, a r eanalysis based 



on th e x-ray detection of lPavlov et al.l ([200a) produced a viable candidate (IMignani et al, 



20081 ). Extrapolating the VLBI position and proper motion back to the date of the opti- 
cal data (July /August 2000) provides further confirmation of the optical detection (VLBI 
position: RA 01:08:08.3 11(1), dec -14 : 31:49 12(1); optical position: RA 01:08:08.301(26), 
dec —14:31:49.15(27)). IMignani et al.l (j2008l ) show that the measured optical emission is 
compatible with purely thermal emission from the bulk of the neutron star, with a surface 
temperature of (7 — 10) x 10^ K, assuming a distance of 130 pc and a neutron star radius 
of 13 km. They note, however, that a determination of the actual neutron star surface 
temperature is dependent on an accurate distance measurement. Taking the lower limit 
of 7 X 10^ K from above, and combining with the 95% confidence lower limit for distance 
to PSR J0108-1431 (143 pc), a simple re-scaling shows that these VLBI measurements 
imply that the temperature of PSR J0108-1431 must be greater than 8.5 x 10^ K. At the 
best-fit distance of 240 pc, the implied temperature is in the range (2.4 — 3.4) x 10^ K. 

Since the temperature of a neutron star with the age of PSR J0108-1431 (Shklovskii- 
corrected characteristic age 200 x 10^ yr) should be less than 3 x 10^ K if it had cooled 
accor ding to standard models without any internal reheating (see Figure 3 of 



Schaab et al 



19991 ). these results imply that some reheating of isolated, slowly rotating neutron stars 



must occur. While PSR J0108-1431 could be younger than its characteristic age, assuming 
standard cooling to a temperature at or higher than the 95% confidence lower limit (8.5 x 
10^ K) implies an age less than 10 x 10^ yr, which is extremely unlikely. More accurate 
temperature measurements for a number of pulsars would be required to differentiate the 
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I999I . for a review of reheating models) but 



possible reheating models (see ISchaab et al 
for PSR J0108-1431 most reheating processes do not have sufficient efficiency. Alternate 
models based on neutron star interiors composed of quark matter s ugges t a slower cooling 
rate and higher temperatures for old neutron stars (jAlford et alJ . l2005l ). but the surface 
temperature implied by these PSR J0108-1431 observations is challe nging to explain even 
in that framework. Alternatively, as noted in iMignani et al.l (j2008l ). a component of the 
optical emission may be non-thermal, although in this case the efficiency of the non- 
thermal process would be higher than in other optically-detected neutron stars such as 



PSR B0950+08 ( Zharikov et al 



200J). 



The predicted velocity of 194 km s~^ is unremarkable, and the pulsar's high Galac- 
tic latitude of —76.8° makes is difficult to place good constraints on the true vertical 
component of velocity. This renders a comparison of characteristic to kinematic age im- 
possible for PSR J0108-1431, but as the pulsar would likely have already completed more 
than one oscillation through the Galactic disk, such a comparison is not useful in any 
event. The proper motion is compatible with, but r nuch more significan t, than the value 



of /^Q = 92 lb 44, fis = 176 ± 70 mas yr ^ derived by 



Pavlov et al 



([2003). 



6.2.2 PSR J0630-2834 



PSR J0630-2834 is a middle- aged pulsar which has been obs erved extensively as part 
of scintillation studies (see e.g. ICorded . Il986l : iBhat et al.l . Il999l ), and would b e undi stin- 



guished if not for its remarkable apparent luminosity in x-rays. 



Becker et al 



((20051) ob- 
served PSR J0630-2834 with XMM-Newton and found an x-ray luminosity of 8.4 x 10^° erg 
s~^, based on the NE2001 distance of 1.45 kpc. This implied that the pulsar was convert- 
ing 16% of its spin-down luminosity into x-rays - an order of magnitude more than any 
other old or middle-aged pulsar. They suggest that the most likely explanation is that 
the distance to the pulsar is over-estimated, a theory which has been proven correct by 
these VLBI observations. 

The measured and fitted positions of PSR J0630-2834 are shown in Figure 16.61 The 
measured distance of 332^^q pc, as shown in Table 16. H means that the actual x-ray 
conversion efficiency is a much less surprising 0.8%. This is within the la error bars of 
the efficiency derived for PSR J0108-1431 in Section 16.2.11 above. The best-fit distance 
of 332 pc is less than a quarter of the DM-derived distance in the NE2001 model, and a 
factor of seven smaller than the distance predicted by the TC93 model (2150 pc). This 
large discrepancy is discussed further in Section [6.31 

The proper motion of PSR J0630-2834 has previously been measured using the VLA 
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Figure 6.6 Motion of PSR J0630-2834, with measured positions overlaid on the best fit. 
From top: Motion in dechnation vs right ascension; right ascension vs time with proper 
motion subtracted; and dechnation vs time with proper motion subtracted. 
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by lBrisken et alJ (l2003l ). and the measured values from this work (//q, = —46.3 it 0.99 mas 
yr~^, fis = 21.26 it 0.52 mas yr^^) are consistent with, but marginahy more precise, than 
the VLA results (/iQ, = —44.6 it 0.9 mas yr~^, fj.^ = 19.5 it 2.2 mas yr~^). The resultant 
Shklovskii correction is less than 0.04% of the observed period derivative (7.12 x 10~^^). 
The improved proper motion accuracy makes a mar ginal impact on the rn easurement of 

(120051 ) , reducing the 



Johnston et al 



the alignment of the rotation and velocity vectors by 
uncertainty in velocity position angle from 3° to 1°. This halves the error in the position 
angle difference measurement and confirms that the velocity and rotation vectors of PSR 
J0 630-2834 appear t o be very highly aligned. 

jwm measured the scintillation velocity of PSR J0630-2834 to b e 170 ± 



Bhat et al 



15km s ^, considerably higher than the earlier measuremen t of 60 =fc 13km s ^ by 
(|l986l ). However, the estimated distance used by 



Cordes 



Bhat et al 



(1199911 wa s the T C93 distance 



(2150 pc), which considerably biases the estimated velocity. ICordesI (Il986l ) uses a pre- 



TC93 d istance estimate o f 124 pc. When re calculated using the VLBI distance of 332 



pc, the 



Bhat et al 



(|l999l ) and ICorded ( 19861 ) speeds are revised to 67 ± 6 km s ^ and 



31 lb 7 km s ^ respectively. The transverse velocit y of 80"*"!^ 



cm s ^ obtained from these 



Bhat et al 



VLBI measurements is thu s consis t ent w ith the 
inconsistent with the earlier ICordea ()1986l ) measurement. 



(jl999l ) measurements, but 



The characteristic age of PSR J0630-2834 is 2.77 x 10^ years, while its kinematic age 
can be calculated from its Galactic latitude b = —16.758° and proper motion in Galactic 
latitude fi^ = —34.8 mas yr~^ as 1.73 x 10^ years. At the best-fit distance of 332 pc, the 
Galactic height of PSR J0630-2834 is —96 pc, and thus the characteristic and kinematic 
ages are consistent, and imply that the progenitor for PSR J0630-2834 resided within 60 
pc of the plane. 



6.2.3 PSR J1559-4438 

PSR J1559~4438 is a relatively unremarkable, middle-aged, isolated pulsar, which was 
included in this observing program as the primary "technique check" source. Being rela- 
tively bright and a strong scintillator, it provided an excellent opportunity to check the 
data reduction tools developed in this thesis in a relatively high SNR regime. However, 
as it has been well studied in the past due to its relative proximity to the Solar System 
and brightness, a number of interesting comparisons to previous studies can be made. 
Although plots of the fitted and measured positions of PSR J 1559-4438 have been shown 
previously in Chapter [Sj they are shown here again in Figure [6771 for completeness. 

The measured distance of 26OOI45Q consistent with the NE2001 prediction (2350 
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Figure 6.7 Motion of PSR J1559-4438, with measured positions overlaid on the best fit. 
From top: Motion in dechnation vs right ascension; right ascension vs time with proper 
motion subtracted; and dechnation vs time with proper motion subtracted. 
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pc iCordes and Lazid . l2002l ) , which differed considerably from the earlier TC93 distance 
estimate of 1580 pc. It is also consist ent with the lower dist ance estimate of 2.0 it 0.5 

(|l995l ). The measured values of 



Koribalski et al 



kpc made using HI line absorptions by 
proper motion (/Iq, = 1.52 it 0. 14 mas yr~^, = 13. 15 it 0.05 mas yr~^) are consistent 

(jl997l ). who measured /Xq, = 1 it 6 mas 



Fomalont et al 



with the VLA observations of 
yr~^, = 14 lb 11 mas yr~^. Neglecting acceleration from the Galactic potential, the 
kinematic age of the pulsar can be estimated from its Galactic latitude b = 6.367° and 
proper motion perpendicular to the Galactic plane fib = 8. 93 mas yr~^ as 2.57 =fc 1.02 
Myr, assuming a birth location within 100 pc of the plane (jFaucher-Giguere and Kaspi . 



20061 ) . Under the standard assumption of a braking in dex of 3, the observe d period 
P = 257ms and period derivative P = 1.01916 x 10^^^ ( Siegman et al. . Il993 ) imply a 
birth pe riod Pn between 35 and 151 ms - longer than is often assumed for normal pulsars 

3 ut sim ilar to the calculated value of Pq = 139.6 ms for 



(see e.g. iMigliazzo et al 



2002 L 



PSR J0538+2817 (IKramer et al 



20031 ) . The correction to the observed P due to the 



Shklovskii effect is only 0.03% of the observed value. 

The transverse velocity measured for PSR J1559-4438 (164 km s^^, 95% confidence 
upper l imit of 287 km s ~^) is inconsistent with the 400 km s~^ estimated from scintilla- 
tion by iJohnston et al.l (|l998l ). who assumed a 2 kpc distance to PSR J1559~4438 and a 
scattering screen midway to the pulsar. If this discrepancy is interpreted as an error in 
the scattering screen location, then the scattering screen must reside considerably closer 
to the pulsar than to the solar system. Alternatively, the distribution of turbulence in the 
scattering disk may be anisotropic, as appears to be the case for PSR J0737~3039A/B. 

With an accurate proper motion now calculated, the position angle of the proper 
motion for PSR J1559-4438 can be compared to the position angle of the emission po- 
larisation, which tests the alignment of the rotation and velocity vector, as described by 
(|2007l ). If the pulsar emits predominantly parallel or perpendicular to the 



Johnston et al 



magnetic field lines, then the angle bet ween the velocity and p olarisation position angles 

([20071) found plausible alignment 



Johnston et al 



is expected to be 0° or 90° respectively, 
in 7/14 pulsars of similar ages to PSR J1559-4438. From Table \67l\ it is easy to calculate 
the velocity position angle as PA„ = 6.6 = b 0.6°. The polarisation position angle for PSR 
J1559-4438 given in I Johnston et al.l ((20071) is 71 ± 3, and so there does not appear to be a 
strong case for alignment of the velocity and rotation axes for PSR J1559-4438. However , 
the polarisation profile of PSR J1559-4438 (as shown in Figure 5 of iJohnston et al.ll2007l ) 
is complicated, and it is possible that the original determination of the magnetic field 
orientation from the polarisation position angle was incorrect. Thus, the alignment of the 
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rotation and velocity vectors for PSR J1559-4438 cannot be excluded, and more detailed 
polarisation studies are required. 



6.2.4 PSR J2048-1616 

PSR J2048-1616 was the second "technique check" source included in the observing pro- 
gram, after PSR J1559-4438. It is less bright than PSR J1559-4438, but was predicted 
to be somewhat closer (640 pc in the TC93 model, 560 pc in the NE2001 model). As 
shown in Table 16.11 whilst a measurement of parallax was made for PSR J2048-1616, it 
was not significant (1.9cr). Nevertheless, the best-fit distance of 580 pc is consistent with 
the DM-based distance estimates, and an accurate measurement of proper motion was 
made. Figure [6^ shows the fitted and measured positions of PSR J2048-1616. 

It is apparent from Figure 16.81 that the large position error of the third epoch (MJD 
54182) is the primary reason for the poorly constrained fit. Equipment failure at the Mo- 
pra telescope and limited time on the Tidbinbilla telescope during this epoch reduced the 
number of telescopes on-source during observations of PSR J2048-1616 to three, and thus 
the large errors are unsurprising. The covariance between proper motion and parallax 
was the largest for PSR J2048-1616, and hence adding even a single future position mea- 
surement at the appropriate parallax extrema would almost certainly allow a significantly 
improved measurement of parallax. 

The measured proper motion of Ha = 114.24 it 0.52 mas yr~^, /i^ = —4.03 it 0.24mas 
yr~^ is consistent wi th but considerably more precise than the VLA measurement of 



Fomalont et al. 



(jl997l ). which was /Xq, = 117 ± 5 mas yr~^, = — 5 ± 5 mas yr"-*^. This 
allows an improved measurement of the transverse velocity vector position angle of 92.0 it 
0. 1 5° . Since the error in t he alignment of polarisation and velocity position angle calculated 



bv I Johnston et al.l (|2007l ) was dominated by the error in polarisation position angle (5°), 
this improvement is not currently useful, although future improvements in polarisation 
measurements will allow a more precise measurement of the position angle difference and 
hence the alignment of the rotation and velocity vectors. The Shklovskii component of P 
for PSR J2048-1616 can be calculated from the proper motion and best-fit distance to be 



3.5 X 10"^^, or 



Bhat et al. 



ess than 0.4% of the observed value of P = 1.0958 x 10"^^. 



(|l999l ) measured the scintillation velocity of PSR J2048-1616 to be 501 ± 
29 km s~^ assuming the TC93 distance of 640 pc, which cannot currently be ruled out 
due to the distance uncertainty. However, if PSR J2048~1616 does indeed reside at the 
best-fit distance of 580 pc (which is consistent with the DM estimates) then, as with PSR 
J1559-4438, the scintillation velocity has been overestimated by a factor of ~ 2, and as 
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Figure 6.8 Motion of PSR J2048-1616, with measured positions overlaid on the best fit. 
Prom top: Motion in dechnation vs right ascension; right ascension vs time with proper 
motion subtracted; and declination vs time with proper motion subtracted. 
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with PSR J1559-4438, this would support the presence of a scattering screen closer to the 
pulsar than to the solar system. 

The characteristic age of PSR J2048-1616 is 2.84 x 10^ yr"^, while the kinematic age 
can be calculated from the Galactic latitude of —33. 077° and proper motion in Galactic 
latitude (-107 mas yr~^) as 1.113 x 10^ yr. As with PSR J1559-4438, assuming a braking 
index of 3 and taking the observed period P = 1961.6 ms and period derivative P = 
1.0958 X 10^^^ allows the calculation of birth period Pq = 1192 it 243 ms, assuming the 
pulsar progenitor resided within 100 pc of the plane. 

Since the distance to PSR J2048-1616 is still quite uncertain from these measurements, 
the most conservative estimates for Pq and/or birth height can be obtained by assuming 
the pulsar distance is the minimum allowed value. At the 95% confidence lower limit of 
distance (283 pc), the Galactic height of PSR J2048-1616 is 155 pc, and the lower limit 
upon Pqj assuming a birth event 100 pc above the plane, becomes 695 ms. Alternately, if 
the pulsar was born with a negligible spin period, the Galactic height of its progenitor 
must have been greater than 240 pc (at 95% confidence). 

While, as noted in Chapter [21 considerable debate still exists as to the true distribution 
of initial pulsar spin periods, values greater than half a second are difficult to obtain with 
standard pulsar birth models. The second alternative, a progenitor residing considerably 
further above the Galactic plane than usual, is possible but unlikely. The remaining 
alternative to reconcile the characteristic and kinematic ages of PSR J2048-1616 is that its 
braking index is greater than the canonical value of 3.0. Substituting the 95% confidence 
upper bound on kinematic age (1.83 x 10^ yr, assuming the minimum distance to PSR 
J2048-1616 and a progenitor 100 pc above the plane) into Equation 12.31 and assuming 



Yueetal 



zero i nitial spin period yields a minimum braking index of 4.1. As noted by 
()2007l ). however, none of the six reliable measurements of braking index made to date 
have been greater than 3.0, which would make such a high value surprising. Clearly, PSR 
J2048-1616 is somewhat unusual in terms of generally assumed pulsar birth and evolution 
parameters, and and further improvements to the parallax distance to PSR J2048-1616 
would be extremely valuable in constraining the birth parameters and spin evolution of 
this pulsar. 
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6.2.5 PSR J2144-3933 



PSR J2144-3933 was discovered in the Parkes Southern Pulsar Surve y ([Manchester 



19961 ). and was initially misidentified with a period of 2.83 seconds. lYoung et al. 



et al 



1999! ^ 



reported that the true period is in fact 8.5 seconds, making PSR J2144-3933 the longest 
period radio pulsar k nown. PSR J2144-3933 lies below the traditionally assumed pulsar 
"death line" (see e.g. 



Chen and Ruderman 



I993I ). and hence its true luminosity is ex- 



tremely important for models of pulsar emission and evolution. The discovery of PSR 
J2144-3933 prompted alternative models of pulsar emission in which long period p ulsars 
such as PSR J2144-3933 could remain luminous in the radio (e.g. lZhang et al.l . 120001 ) . PSR 
J2144-3933 possesses a steeper spectral index than average (—2.4), and while several other 
pulsars appear less luminous at 400 MHz, PSR J2144-3933 is the least luminous pulsar 
observed at 1400 MHz (apparent luminosity 24 /iJy kpc^ at assumed distance of 180 pc; 
average 1400 MHz flux density of 0.75 mJy calculated from archival Parkes observations). 

The fitted and measured positions of PSR J2144-3933 are shown in Figure 16.91 As 
shown in Table a highly significant parallax was detected, corresponding to a distance 
of 165^^4 pc. Given the generally assumed errors on DM distances, this is consistent 
with the TC93 value (180 pc), but not the NE2001 value (264 pc). This confirms that 
the apparent radio luminosity of PSR J2144-3933 is extremely low - 15% lower than the 
previously assumed value of 24 /iJy kpc^. 

The proper motion measurement allows for the first time a calculation of the Shklovskii 
correction for the period derivative of PSR J2 144-3933 - the Shklovskii contribution to 
P is 9.4 X 10^^'^, or approximately 19% of the observed P of 4.96 x 10~^^. The true 
spin-down luminosity of PSR J2144-3933 is thus further reduced from the assumed value 
of 3.2 X 10^^ erg s~^ to 2.6 x 10^^ erg s~^, the smallest known spin-down lum inosity of any 



pulsar. PSR J0343-3000, discovered in the Parkes High-Latitude survey bv lBurgav et al 



(|2006l ). has the next lowest spin-down luminosity, which is a factor of 5 higher than the 
revised value for PSR J2 144-3933. In addition, this revision to the P value for PSR 
J2144-3933 places it even further past the assumed pulsar "death line". 

The correlation between observed pulsar radio luminositie s (L) and their sp in-down 
luminosity {E) has been the subject of considerable debate. Lvne et all (1998) argued 



for a model in which t he pulsar luminosity was ind e 



pendent of the other known 



20061 ) and 



p hysica l 



Malov and Malov 



(|200a). 



pulsar parameters, but iFaucher-Giguere and Kasp] ( 
amongst others, have shown evidence for a correlation with the pulsar's spin-down lumi- 
nosity, with a dependence ranging from L oc E ^^'^ to L oc E^/"^ . More general mod els with 



a dependence upon P and P are considered by iFaucher-Giguere and Kaspil (|2006l ) but are 
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Figure 6.9 Motion of PSR J2144-3933, with measured positions overlaid on the best fit. 
From top: Motion in dechnation vs right ascension; right ascension vs time with proper 
motion subtracted; and dechnation vs time with proper motion subtracted. 
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not investigated here. Since the radio luminosity of ordinary pulsars is believed to be 
powered by their rotational energy, it is logical to assume that the radio and spin-down 
luminosities are related, but as already noted, complications such as differing beaming 
geometries make it plausible that this relationship could be difficult to identify observa- 
tionally. Additionally, since the radio emission of pulsars makes up only a small component 
of their total emission, there is a potential for a large scatter in any relationship between 
radio luminosity and spin-down luminosity, and pulsar spectral indices can vary consid- 
erably, meaning measurements based on a single frequency introduce a further source of 
uncertainty. Certainly, spin-down luminosity and observed 1400 MHz radio luminosity 
appears to be uncorrelated when the entire population of pulsars is considered, as shown 
in the left panel of Figure I6.10[ 

However, it is interesting to note that PSR J2144-3933 has both the lowest spin- 
down luminosity a r id the lowest apparent 1400 MHz radio luminosity of all known pulsars. 



Malov and Malovl (|2006l ) show that the radio efficiency of old pulsars (P > 0.1 s) appears 
to increase with period, with the conversion efficiency for 2 second pulsars on average being 
an order of magnitude greater than for pulsars with period 0.2 seconds. PSR J2144-3933 
is the longes t perio d pulsar known and challenges standard pulsar emission models (e.g 



Zhang et al.l . l200d ). so it is reasonable to investigate the possibility that pulsars with 



similar characteristics might show a stronger correlation between s pin-down luminosity 



and r adio luminosity. The vacuum gap model of pulsar emission (e.g. iBhattacharva et al 



19921 ) predicts a radio death line given by: 

B/P^ < 1.7 X 10"G s"2 (6.6) 



The right panel of Figure 16.101 plots apparent radio luminosity against spin-down 
luminosity for pulsars that lie near or below this radio death line. Since only 21 pulsars 
lie below this death line, and 1400 MHz luminosities are available for only 14 of these, the 
B / cut-off was chosen to be twice the vacuum gap death line value (3.4 x 10^^ G s"^). 
This resulted in a sample of 65 pulsars. 

As shown in Figure 16.101 the correlation between spin-down luminosity and apparent 
1400 MHz radio luminosity is stronger for pulsars that lie near the traditional pulsar death 
valley, towards the lower right corner of the P-P diagram. The line of best fit shown in 
the right panel of Figure 16.101 was obtained using standard linear regression, and is given 
by the equation log^g L = 0.671og^Q E — 19.57. The 95% confidence interval for the power- 
law exponent is [0.24, 1.1], and the value of 0.136 is inconsistent with no correlation at 
greater than 99% confidence for this sample size of 65 pulsars. These results support the 
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Figure 6.10 Correlation between apparent radio luminosity and spin-down luminosity for 
all pulsars (left panel), and pulsars with periods near or below the traditional death line 
[B/P'^ < 3.4 X lO^^G/s; right panel). Whilst neither correlation is particularly strong, 
the value is higher for the "death valley" pulsars in the right panel (0.136) compared 
to the whole population shown in the right panel (0.033). 



L oc E^'^ model presented by 



Malov and Malov 



(j2006l ). but models with a considerably 



steeper or shallower exponent cannot be excluded. 

As already noted, the unknown effect of the beaming geometry for each pulsar is bound 
to impose a considerable scatter upon the apparent radio luminosities. Furthermore, all 
spin-down luminosities are uncorrected for the Shklovskii effect and, as discussed further 
in Section 16.3.21 errors in pulsar distances could be corrupting the radio luminosities 
and further obscuring a more obvious correlation. The results presented here should 
not be interpreted as definitive evidence of a correlation between integrated pulsar radio 
luminosities and spin-down luminosities - rather, they suggest that such a correlation is 
more likely for pulsars near the death valley, and that more detailed m odeling could help 



Malov and Malov 



confir m or deny the existence of a correlation such as that proposed by 
(|200g). Such modeling would necessarily include obtaining accurate distances and spectral 
indices and modeling beam geometries to estimate the integrated radio luminosity, and 
obtaining proper motions to calculate Shklovskii correctio ns to spin-down lu minosity. 

As noted in the discovery paper of PSR J2144-3933 bv I Young et al.l (jl999l ). the low lu- 
minosity and narrow beaming fraction (~ 0.01) of PSR J2144-3933 imply that many such 
objects exist in the Galaxy. The distance measurement provided by this VLBI astrometric 
program has confirmed and slightly strengthened this hypothesis, since the pulsar is even 
closer and fainter than previously thought. Assuming the steep spectral index of PSR 
J2144-3933 is typical for similar long-period, low-luminosity objects, and that the Galac- 
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tic population of such objects c an be estimated fro m PSR J2144-3933 alone (which yields 
a population of 100,000 objects: lYoung et alJll999l ). many such objects would be expected 
to be detected by LOFARo. LOFAR, which is a low-frequency aperture array telescope 
presently being constructed in the Netherlands, will search for pulsars at frequencies of 
approximately 120 MHz. 

While the spectrum of PSR J2144-3933-like objects may exhibit a turn-over before 
this point, assuming the measured spectral index of —2.4 continues to 120 MHz yields 
a 120 MHz radio luminosity of 8 mJy kpc^. While the scope of LOFAR has been ad- 
justed several times recently, in pulsar survey mode it is expected to have a point source 
sensitivity in the 120 MHz band of ^ 0.2mJ}|^. Thus, LOFAR would be expected to 
detect objects like PSR J2144~3933 up to distances of several kpc, beyond which distance 
interstellar scattering is likely to limit detections more than sensitivity. Thus, potentially 
thousands of PSR J2144-3933-like objects could be expected to be detected with LOFAR. 
Low-frequency observations of PSR J2 144-3933 would resolve the question of a potential 
spectral turnover, and help refine estimates of LOFAR detection rates. 

W hile this is the first measurement of proper motion for PSR J2144-3933, iJohnston et al. 



( 19981 ) have measured a scintillation velocity of 48 km s~^ assuming a distance of 330 pc, 
considerably lower than the VLBI-derived value of 130^^2 Since the scintillation 

speed was underestimated, despite the fact that the distance to the pulsar was overes- 
timated, a possible explanation is that the scattering screen resides closer to the solar 
system than the pulsar. Alternatively, the scintillation measurement could have been af- 
fected by refractive scintillation, or the thin-screen, isotropic turbulence model may not 
be applicable in this case. 

No detections of PSR J2144-3933 in optical or x-ray wavebands have been published. 
With a Shklovskii-corrected characteristic age of 336 x 10® years, it is even older than 
PSR J0108-1431, and might be expected to possess similarly high conversion efficiency 
from spin-down luminosity to optical and x-ray luminosity. However, even though it is 
considerably closer than PSR J0108-1431, the fact that its spin-down luminosity is a factor 
of almost 50 smaller means that it would be challenging to detect in wavebands other than 
radio. A simple scaling of the optical observations of PSR J0108-1431 to account for the 
different distance and spin-down luminosity of PSR J2144-3933 reveals that the U band 
optical magnitude would be 27.7. which would be detectable with large ground-based 
telescopes, although challenging. The x-ray luminosity, assuming the 1.7% conversion 



http://www.lofar.org/ 

''assuming a five minute integration, 24 Mffz bandwidth, and a 10% pulsar duty cycle - see 
http://www.lofar.Org/p/astronomy_spec.htm 
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efficienc y seen in PSR J0108-1431, would be 5.5 x 10^^ erg s"^ Scaling the count rate 
seen by 



Pavlov et al 



((20081) on PSR J0108-1431 by the relative distance and luminosity 
of PSR J2144-3933 compared to PSR J0108-1431 yields the extremely low photon count 
rate of 0.08 counts ks^^, and thus x-ray detection of PSR J2144-3933 is unlikely without 
an extremely long exposure. 



6.3 Analysis of distance and velocity models 



Galactic electron distribution models are extremely important for all fields of pulsar re- 
search, since their prediction of pulsar distances based on DM has the potential to bias 
estimates of luminosity, velocity, and various timing terms for the vast majority of pulsars 
without an independent distance constraint. Hence, continual improvement of these mod- 
els, and characterisation of errors in the existing models (and the impact of these errors) 
is a crucial task. 

Early attempts to construct si mple Galactic electron models fo r the purpose of estimat 
ing pulsar distances were made by iManchester and Taylor! (jl98ll ) and 
but th ese were superseded by the more comprehensive T C93 model of 
( I993I ). While the much more complex NE2001 model of 



Lvne et al.l (Il98,5l) 



Tav 



Cordes and Lazio ( 



or and Cordes 



2002 



IS now 



available, distance estimates based on the TC93 model are still in common usage and both 
models are considered below. 

The six new (and one refined) pulsar parallaxes measured in this thesis make a sub- 
stantial addition to the 29 published (18 VLBI and 11 timing) pulsar parallaxejfl and hence 
a review of the accuracies of the TC93 and NE2001 models is timely. It is appropriate to 
note at the outset of this analysis that large distance errors may be over-represented in 
existing astrometric results. One reason is the selection effect of anomalous pulsars being 
chosen for study, which is certainly the case for the PSR J0630-2834 system observed 
during this thesis. Another potential factor is that most astrometric distance determina- 
tions to date have been for relatively nearby pulsars, where there is less opportunity for 
underdensities and overdensities in the ISM along the line of sight to cancel, and hence it 
is plausible that distance models are on average more reliable for more distant pulsars. 

Figure 16.111 plots the TC93 distance for all 35 pulsars against the parallax distance, 
along with a line of best fit. Figure [6.121 repeats this plot for NE2001 distances. The best 
fit line was generated by averaging the linear regression results obtained from a 10,000 
trial Monte Carlo simulation, where "actual" pulsar distances were calculated for each trial 



* Only the most accurate parallax measurement has been considered where multiple measurements 
exist, and measurements less significant than 1.5ct have been excluded 
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Figure 6.11 Parallax distance versus TC93 distance for all pulsars with a published paral- 
lax. The line of best fit was estimated from a Monte Carlo simulation of the true pulsar 
distances based on parallax measurements and errors, averaging the results obtained from 
standard linear regression. The TC 93 model is particularly uncertain for nearby pulsars, 
with underestimates and overestimates by factors of up to six. 



based on the observed parallax and parallax error. Throughout this section, errors will be 
referred to in decibels (dB), which is the most convenient representation for measurements 
where large underestimates or overestimates are common. The error in dB is defined as 
10 logxo p™r'an'ax'^distance ' hence positive values represent overestimates, and negative 
values represent underestimates. 

Figures 16.111 and 16.121 show that the NE2001 model has considerably improved upon 
the TC93 model, but significant errors are still made for individual pulsars. This can 
be highlighted more clearly by binning the errors, as shown in Figure 16.131 for the TC93 
model and Figure [6. 141 for the NE2001 model. Figure [6.131 shows a clear systematic offset 
in the median error of the TC93 model towards underestimated distances, but the largest 
errors are seen when distances are overestimated. Figure 16.141 shows that systematic bias 
has been removed by the NE2001 model, but the distribution of errors still cannot be well 
approximated by a single Gaussian, with long wings towards large errors. For both the 
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Figure 6.12 Parallax distance versus NE2001 distance for all pulsars with a published 
parallax. As with Figure 16. IH the line of best fit was estimated using linear regression 
in a Monte Carlo simulation to account for parallax measurement errors. The RMS error 
is substantially reduced compared to the TC93 predictions, but as some pulsar distances 
had been measured prior to 2001 and were used to constrain the NE2001 model, part of 
the observed improvement must be attributed to this knowledge. 
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Error (dB) 



Figure 6.13 Histogram of TC93 errors for pulsars with measured parallaxes, with errors 
binned in 1 dB increments. Underestimates of distance are more common in this model, 
but the largest errors are made when the distance is overestimated. The standard deviation 
of the errors is approximately 3 dB. 

TC93 and NE2001 models, the largest errors are seen when the distance is overestimated. 
The effect of these high-error "outliers" on population studies of neutron stars is considered 
further in Section 16.3. 2[ 

6.3.1 Analysis of newly— measured pulsars 

Table [6^31 shows the measured distances and velocities for the pulsars in this survey, along 
with previous estimates of their distance from the TC93 and NE2001 electron models, and 
(where available) estimates of their speed from scintillation studies. PSR J0437-4715 has 
been excluded from the analysis of DM distance reliability below, since it contributed to 
the NE2001 model. Since no parallax was detected for PSR J2145-0750, the VLBI velocity 
is shown assuming the NE2001 distance of 566 pc, but given the consequent uncertainty 
it is not included in the error analysis of scintillation velocity estimates. Naturally, PSR 
J2145~0750 cannot contribute to the estimates of DM distance reliability. 
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Figure 6.14 Histogram of NE2001 errors for pulsars with measured parallaxes, with errors 
binned in 1 dB increments. Unlike the TC93 model, there are no discernable systematic 
errors, but the distribution is still clearly non-Gaussian, with a long tail of errors up to 
6 dB. The occurrence of errors with a magnitude greater than 3 dB has been somewhat 
reduced from the TC93 model (five, compared to eight for TC93). The standard deviation 
of the errors is approximately 2 dB, compared to the 3 dB seen for TC93. 
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Table 6.3. Comparison of distance and velocity to non VLBI-derived estimates 



Pulsar 


VLBI distance 


VLBI velocity 


TC93 dist 


ancc 


NE2001 distance 


Scintillation speed 




(pc) 


{km s^l) 


(pc) 




(pc) 


(km s^l) 


PSR J0108-1431 
PSR J0437-4715 
PSR J0630-2834 
PSR J0737-3039A/B 


2401-" 

156.31JI 
332^^^ 

1150+220 


19411- 
104.7+;-° 

80+ ;5 




130 
140 
2150 
570 


184 
189 
1444 
483 


170^ 231^ 
60 ± 13? 170 ± 15° 
140.9 ± 6.2? 66 ± is'' 


PSR J1559-4438 
PSR J2048-1616 
PSR J2144-3933 




1301- 




1580 
640 
180 


2384 
562 
264 


400* 
501 ± 29° 
48* 


PSR J2145-0750 




46t« 




500 


566 


31 ± 25p 51^ 113^ 



* IJohnston et al.l 119981^ using TC93 distances, except for PSR J1559-4438 where 2 kpc (HI lower limit) was used and PSR 
J2144— 3933 where 330 pc (unknown source) was used 

^Gothoskar and Guptal 120001 ). using TC93 distances 

' ^Cordeii 11989) . using a distance of 1240 pc 

' iBhat et al.l 119991 ). using TC93 distances 

' ^Ransom et ahl ]2004l ) . using a distance of 600 pc 

es et al.l t2005l) . using a distance of 500 pc 

nNicastro and Johnstonl 119951 ). using the TC93 distance 



Table ES] shows that for this admittedly limited sample size, the NE2001 model is only 
slightly improved from the TC93 model, with a mean distance error of 2.1 dB compared 
to 2.5 dB. This is smaller than is seen for the entire population, where the NE2001 model 
shows a 1 dB improvement over the TC93 model. This may be partly due to the presence 
of PSR J0630-2834, which dominates the error budget for both the TC93 and NE2001 
models, as the models overestimated the distance to this pulsar by a factor of 6.5 and 4.3 
respectively. Clearly, for individual pulsars, DM distances alone cannot be relied upon to 
provide luminosity, velocity or space density estimates. 

Table [6^31 shows that scintillation speed estimates are also generally unreliable for any 
given system ~ of the nine estimates available, only two are consistent with the VLBI- 
derived velocities, and one of those (PSR J2048-1616) is likely to also be inconsistent, but 
cannot presently be ruled out due to the low significance of the parallax detection. In 
some cases (such as PSR J0630-2834) the incorrect distance estimate to the pulsar was 
the main reason for the scintillation speed estimate errors, but for most cases the DM 
distance error did not significantly contribute to the scintillation speed error. The mean 
absolute error is 3.8 dB, with the largest error (on PSR J0737-3039A/B) being a factor of 
7. Amongst this limited sample, scintillation speeds have consistently overestimated the 
true tra nsverse speed, on average by a factor of 1.75. This is in contrast to previous authors 
such as iGuptal ()l995l ). who found that scintillation speeds tended to underestimate the 



proper motion speed, at least for pulsars with a large Galactic height. These findings are 
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little changed when the scintillation speed estimates are corrected by the VLBI distance 
measurements - the mean overestimate of speed becomes a factor of 1.5, and the mean 
absolute error is 4.1 dB. 

Most previous studies of scintillation speed have assumed a single, thin, turbulent 
scattering screen midway between the pulsar and the Solar System. The turbulence within 
the screen is generally assumed to be isotropic. Recent results from scintillation studies of 
PSR J0737-3039A/B (W. Coles et al., in preparation) and PSR B0834+06 (W. Brisken 
et al., in preparation) suggest that anisotropy of ISM turbulence is the norm, rather than 
the exception. When coupled with the uncertainty of the scattering screen location, this 
suggests that errors in scintillation speeds have typically been underestimated in the past. 
Whilst the scintillation speed method remains a valid approach for obtaining population 
velocity measurements (subject to the same caveats highlighted in Section I6.3.2p . the 
natural conclusion is that too many uncertainties exist for the method to be useful for 
individual objects, at least when no information on the location and structure of the 
scattering material is available. 

6.3.2 Impact of DM distance errors 

It has already been shown that although there is no evidence for systematic bias in the 
NE2001 distance model, errors exceeding 4 dB still exist for individual pulsars. Such errors 
can impact upon population analysis of pulsars, artificially creating tails of extremely high 
or low values of distance-dependent parameters such as luminosity and velocity. This can 
affect estimates of the mean values of these quantities, as well as confusing attempts to 
explain the underlying physics by generating false outliers. As an example, the case of 
neutron star transverse velocities is considered below. 

In order to model the effect of distance errors on transverse velocity estimates, it 
is necessary to estimate the form of the error distribution shown in Figure 16.141 Since 
the small number of samples makes an estimation of the true distribution difficult, the 
error distribution has instead been modeled by binning the errors in increments of 0.1 
dB and smoothing with a 2 dB Hanning window. This is shown overlaid on the original 
binned error distribution in Figure [5.3.21 and is hereafter referred to as the binned error 
model. As already noted in Section 16.31 large errors are likely to be somewhat over- 
represented due to selection effects, and thus this model can safely be considered to be a 
"worst-case" representation of DM distance reliability. Figure 16.3.21 also shows a single 
Gaussian component model for the error distribution, with a standard deviation of 0.66 
dB. In this Gaussian distribution, 75% of errors are smaller than a factor of 1.2, and hence 
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Error (dB) 

Figure 6.15 Histogram of NE2001 errors for pulsars with measured parallaxes, with error 
binned in 1 dB increments. The "binned" error model described in the text is shown as 
a dashed red line, while the single Gaussian "traditional" model which approximates the 
20% errors commonly assumed for DM distances is shown as a yellow dashed line. 



it approximates the 20% errors typically assumed for DM models. This distribution is 
hereafter referred to as the traditional error model. 



Hobbs et al 



(|2005l ) examine a large sample of pulsar proper motions and conclude 
that the distribution of 3D space velocities of young pulsars (characteristic age < 3 Myr) 
is well fit by a Maxwellian distribution with mean 431 km s"^ and one-dimensional RMS 
265 km s^^. This velocity distribution is used as the starting point for the simulations 
that follow. 

A Monte Carlo simulation was performed, creating a synthetic pulsar catalogue of 
1 million pulsars, where each pulsar's actual velocity was drawn from the distribution 
described above. These are referred to as the unperturbed velocities. The velocities of each 
pulsar were then perturbed according to the binned error function and the traditional error 
function described above, and the observed 2D velocity was recorded for each case. The 
observed 2D velocity distribution for the binned error function and the traditional error 
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Figure 6.16 2D velocities for a synthetic pulsar population, as observed with no distance 
error (red), the traditional distance error model (blue) and the binned distance error model 
described in the text (green). The effect of the binned distance error model is to broaden 
the observed velocity distribution, particularly at large velocities. 



function, along with the unperturbed 2D velocity distribution, is plotted in Figure [6.161 



It is immediately apparent that the broadening of the observed velocity distribution is 
considerably greater for the binned error model, compared to the traditional error model, 
which is particularly pronounced at the high end of the velocity distribution. The mean 
unperturbed 2D velocity is 332 km s~^, which is increased by only 1% to 336 km s~^ in 
the traditional error model, but by 18% to 392 km s~^ in the binned error model. This 
in crease exceed s the ± 40 km s"^ errors quoted for the mean 2D velocities of young pulsars 

mm- 



m 



Hobbs et al 



Even more important than the mean effect on pulsar velocities, ho wever, is the effect 
on the high end of the pulsar velocity distribution. iHobbs et al.l (j2005l ) note that distance 
model inaccuracies could be responsible for the sparse tail of very high velocity pulsars, but 
conclude that with the small sample size available, the observed high-velocity pulsars are 
consistent with the tail of continuous velocity distribution. In considering the impact of 
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distance model errors, it is useful to examine the frequency of occurrence of pulsars with 
observed transverse speeds greater than 1000 km s^^ in the synthetic catalogue. These 
pulsars will be described as very high velocity (VHV) pulsars. The synthetic catalogue 
contains only one VHV pulsar in every 1000 in the unperturbed velocity distribution, but 
the occurrence rises to three in 1000 for the traditional error model, and 50 in 1000 for 
the binned error model. 

Thus, distance model errors can add a significant high-velocity tail to the observed 
pulsar velocity distribution, an effect which is significantly enhanced when the distance 
errors "outliers" are accounted for. These results suggest that the mean pulsar 3D space 
velocity may have been previously overestimated, and that the majority of VHV pulsars 
may in fact simply be pulsars with incorrect distance estimates. As already stated, the 
binned error model used here is likely to overestimate the true frequency of large errors, 
and hence the true impact of DAd distance errors is likely to lie somewhere in between the 
binned error model and the traditional error model. Additionall y, there is clear evidence 
that at least some VHV pulsars are present in the Galaxy, with IChatterjee et al.l (|2005l ) 
using VLBA astrometry to show that the transverse velocity of PSR B1508+55 exceeds 
1000 km s~^, and observations of PSR B2224+65 showing scintillation, proper niotion 
and bow shock resu l ts consistent with a velocity close to 1000 km s~^ ( Cordesl . 119861 : 



Harrison et al. 



1993 



Cordes et al 



19931 ). Nonetheless, the binned error model is useful 



in illustrating that the largest impact of occasional, rare distance errors is at the high 
end of the pulsar velocity distribution. A larger, unbiased sample of pulsar distance 
measurements would enable a more accurate estimation of the true distance error function 
for NE2001, which in turn would enable the effects of distance errors to be "deconvolved" 
from the measured pulsar velocity distribution with some confidence. 



6.4 Systematic errors and astrometric accuracy 



Previous astrometric VLBl programs such as those carried out bv lChatteriee et al.l (j2004l ) 
have attempted to estimate the effect of various observing factors on the resultant astro- 
metric accuracy. Two of the most important factors are the brightness of the pulsar, which 
determines the accuracy to which a centroid position can be estimated within the synthe- 
sised beam, and the angular separation of the pulsar from the calibrator (the "calibrator 
throw"), which determines the extent to which atmospheric and ionospheric gradients 
contribute to systematic errors. Broadly, these could be considered as dominating the 
thermal (random) error component, and the systematic (non-random) error component 
to a position measurement. It is appropriate to note at the outset of this analysis that the 
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Figure 6.17 Parallax error plotted against ungated pulsar flux density (mJy, 1400 MHz) for 
all pulsars with a published VLBI parallax. Whilst stronger pulsars typically have smaller 
errors, there is a large scatter. The three faintest pulsars with measured parallaxes to date 
were part of this observing program. 



small number of measured parallaxes and the potential for selection biases (for example, 
some "easier" pulsars have been targeted by earlier surveys with less sensitivity, leading 
to larger errors compared to what might be expected) mean the trends presented here 
should be interpreted with caution. 

Figure 16.171 shows the ungated pulsar flux density (at 1400 MHz) plotted against 
parallax errors for all pulsars with a published VLBI parallax. No attempt was made to 
account for the effect of varying calibrator throw. Where several parallax measurements 
have been made, the most accurate is used. As expected, considerable variation is seen, 
but there is a general trend towards lower parallax errors for brighter pulsars. 

Figure 16.181 shows the parallax errors plotted against calibrator throw for all pulsars 
with a published VLBI parallax. No attempt was made to correct for pulsar flux. As 
with Figure 16.171 the most accurate measurement was taken when multiple measurements 
where available. Where multiple calibrators were used, the shortest calibrator throw was 
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taken (including in-beam calibrators when used). Surprisingly little trend is seen with 
decreasin g calibrator throw, whic h differs from predominantly linea r decrease seen in th e 
results of 



Pradel et al 



(200 



3). 



Chatteriee et alj (|2004l ). and the simulation predictions of 
As previously stated, the sample used here consists of measurements made with many 
different VLBI arrays, at different recording rates and hen ce sensitivities. This approach 
differs from previous studies such as lChatteriee et al.l (|2004l ). who considered single objects 
with different calibrators, or multiple objects with different calibrator throws observed 
with the same array. The in-beam calibrators used also vary considerably in brightness, 
which may limit the attainable precision to less than the expected value from the spatial 
interpolation alone. The variation of pulsar flux amongst the different target objects, 
however, is likely to be the most important factor obscuring the expected correlation 
of parallax error with calibrator throw. This would imply that the majority of pulsar 
astrometry programs to date have been sensitivity-limited, rather than being limited by 
systematic errors. 

Ideally, a joint analysis would be made of pulsar flux and calibrator throw on the par- 
allax error, but the current sample size is too small for such an analysis to be meaningful. 
More tightly controlled and larger samples of parallax measurements will become available 
in future years, enabling a better analysis of the relative contribution of calibrator throw 
to parallax error. 

The pulsars observed in this thesis varied widely in flux and calibrator throw, and so 
different limiting factors are seen for the different targets. As shown in Figure 16.171 the 
three faintest pulsars with measured VLBI parallaxes were observed in this thesis, and the 
results for these three are partially (PSR J2144-3933) or predominantly (PSR J0108-1431, 
PSR J0737-3039A/B) sensitivity-limited. PSR J2145-0750, which was only detected in 
two epochs, was also severely sensitivity limited. Thus, half of the target sample was 
dominated by thermal, rather than systematic errors. For the remaining four pulsars 
(PSR J0437-4715, PSR J0630-2834, PSR J1559-4438 and PSR J2048-1616) systematic 
errors dominated thermal errors and equally weighted visibilities were used to minimise 
the impact of systematic errors on the fitted positions, as discussed in Section [5.5. jpl . 

The systematic error budget is likely to be dominated by residual ionospheric errors, 
since as noted in Section [5.3.21 the relatively low density of GPS receivers in the Southern 
Hemisphere makes the GPS-based ionospheric correction less accurate for the LBA than 
for the VLBA or EVN. As shown in Section 15.5.21 the median position shift due to the 



^PSR J2144-3933, which was moderately sensitivity-limited, used equally weighted visibilities, but 
used natural weighting for imaging to maximise sensitivity, rather than the uniform weighting used for the 
systematic-dominated pulsars 
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Figure 6.18 Parallax error plotted against calibrator throw in degrees, for all pulsars with 
a published VLBI parallax. Little trend is apparent with decreasing calibrator throw for 
both this work and previous surveys, which is surprising given the assumed dependence 
of systematic errors on calibrator throw. 



application of ionospheric correction is approximately 3 mas for PSR J1559-4438. Since 
the accuracy o f GPS -derived TEC maps ranges from 3-10% in the Northern Hemisphere 



( Sekido et al 



20041 ) ■ the accuracy in the Southern Hemisphere is presumably ~ 10% or 



worse, and hence systematic errors due to residual ionospheric effects would be expected 
to be hundreds of //as for pulsars observed at 1600 MHz, comparable to the observed 
systematic error estimates. As discussed in Section I6.1.H the application of ionospheric 
corrections made no significant difference for the 8 GHz observations, and residual tropo- 
spheric errors are expected to dominate in this case. 
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6.4.1 The future of pulsar astrometry 

Over the next decade, the construction of new radiotelescopes and the upgrading of exist- 
ing interferometers to larger bandwidths will lead to a revolution in astrometry at radio 
wavelengths. The most ambitious, and most distant, project is the Square Kilometre Array 
(SKA0), which when completed will combine nanoJansky sensitivity with excellent spa- 
tial resolutiorl^ of ~1-10 mas (in the high frequency band). LOFAR, which has already 
been discussed in Section 16. 2. 5^ should be completed much earlier than the SKA and will 
offer excellent sensitivity at low frequencies, providing astrometric capabilities for faint, 
nearby objects, especially steep-spectrum pulsars. Given its low observing frequencies, 
LOFAR will suffer from relatively low resolution and increased ionospheric effects. For 
the faintest steep-spectrum objects, however, LOFAR will enable astrometry that would 
not otherwise have been possible. 

In addition to these new facilities, large-scale bandwidth upgrades to the VLAp^. 



MERLIN^i. and VLBAo intruments are bringing //Jy sensitivity to instruments with 
baselines of tens, hundreds and thousands of kilometres respectively. Due to its shorter 
baselines, the VLA is unlikely to be able to compete in terms of astrometric accuracy 
generally, but may be useful for high frequency observations of flat -spectrum sources. 

The increased sensitivity of these instruments will allow fainter pulsars to be targeted, 
but will not improve astrometric accuracy for brighter pulsars beyond what can presently 
be obtained unless improved calibration methods can be applied. As shown in Section [6. 1.11 
for PSR J0437-4715, systematic errors already dominate the error budget with existing 
instrumentation for the brightest pulsars, even with the comparatively modest sensitivity 
of the existing LBA, and at frequencies where the influence of the ionosphere are minimal. 

One approach which has already been successfully demonstrated is the use of wide 
bandwidths to measure the quadratically dependent (upon frequency) compone nt of visi- 



bility phase induced by the ionosphere, using self-calibration on the target source (iBrisken et al 



20021 ) . Whilst this technique is likely to be particularly valuable for pulsar astrometry with 
future instruments, it cannot account for tropospheric errors and will still not be feasible 
for the faintest pulsars. Thus, new methods of calibration will be required. 

Ideally, calibration will minimise the error in the estimated path delay and complex 
gain in the direction of the target source at all times, across the entire observing band- 



""^"http:/ /www. skatelescope.org/ 

^^in its current specification, which includes baselines to 3000 km 



^■^http:/ /www.vla.nrao.edu/ 
^^http:/ /www.merlin. ac.uk/ 
^^http:/ /www.vlba.nrao.edu/ 
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width. This removal of the temporal and spatial interpolation which has traditionally been 
required in VLBI phase-referencing experiments is possible using a technique which is com- 
monly known as "field calibration" , which is analogous to the technique of adaptive optics 
in optical astronomy. Whilst the coherence time is much longer at radio wavelengths than 
optical, the density of strong radio sources is much lower than optical sources, meaning 
that this technique could not be applied to the narrow-field, low sensitivity VLBI arrays 
which have been available to date. The advent of wide-bandwidth systems will provide 
the sensitivity required t o make this techn i que p ossible on the planned and upgraded fa- 
cilities mentioned above. iBhatnagar et al.l ( 20081 ) give an overview of the algorithms used 
for field calibration. 

Thus, for astrometric purposes, the availability of more calibrator sources is potentially 
the largest advantage of the improved sensitivity of new and upgraded instruments. Of 
the instruments listed above, the upgraded MERLIN (eMERLIN) is unique in that it is 
a heterogeneous array - all of the other arrays are composed of identical elements. The 
presence of several large telescopes in the eMERLIN array will make the field calibration 
described above more difficult, since the field of view of the array is limited by these large 
telescopes. Thus, it will be difficult for eMERLIN to reach the limiting accuracies that 
should be possible with the upgraded VLBA. 

Since field-based calibration routines have not yet been tested in detail, it is difficult 
to predict the limiting accuracy that can be attained. Ultimately, the ability to model 
the (time-variable) structure of e ach source used as an in- beam calibrator is likely to 
determine the precision achieved. iFomalont and ReidI (|2004l ) estimate that the SKA can 
obtain single-epoch accuracies of several /xas for mJ y targets, meaning sub-//a s parallax 
accuracies should be attainable. However, as noted in Fomalont and Reid 



(j200J), the final 

specification of the SKA instrument, in particular the amount of collecting area devoted 
to long baseline stations, could significantly impact these estimates. The upgraded VLBA, 
unlike the SKA, will still be sensitivity limited for mJy sources - nevertheless, /ias parallax 
accuracies should be attainable, at least for some sources. 

The accuracies for these radio interferometers can be contrasted with the next gener- 
ation of optical instrumentation. The dedicated orbiting astrometric telescope GAIA0, 
scheduled to launch in 2011, will derive accurate positions for over a billion stars. The 
astrometric accuracy will depend on target brightness, but is expected to range from ~7 
/ias for bright sta rs, to ^ 25 Mas for 15th magnitude stars, and sub-mas parallaxes for 20th 
magnitude stars (jLindegren et al.l . 120081 ) . This will allow significant distance determina- 



www.esa.int/science/gaia 
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tion for these stars at distances as great as the Galactic centre. Clearly, GAIA will allow 
the measurement of a much greater number of sources than will be possible with VLBI, 
even in the SKA era. However, the most accurate determination of individual objects is 
almost certain to still be achieved using VLBI. 

Since GAIA will be limited to objects much brighter than pulsars, it will not be 
able to contribute to pulsar astrometry specifically except for biriary p ulsars with an 
optically visible companion, such as PSR B1259-63 (jjohnston et al.l . Il994l ). While many 
pulsar companions are faint white dwarf stars that would be challenging targets for optical 
astrometry, the detection of even a few of the brightest optical pulsar companions would 
enable an improved frame tie between the optical reference frame and the ICRF. Since 
the frame tie between the Solar System frame and the ICRF is already good, thanks 
to comparisons of VLBI and timing astrometry, this would also serve to enhance the 
compatibility of the optical and Solar System frames. 



7 

CONCLUSIONS 



As described in Chapter 1, this thesis aimed to develop the tools necessary to enable 
high-sensitivity pulsar observations with the Australian Long Baseline Array, allowing 
the largest pulsar astrometry program in the Southern Hemisphere to date to be under- 
taken. These tools included a general purpose software correlator (DiFX) and a flexible 
astrometric data reduction pipeline which implements several novel techniques for pulsar 
VLBI astrometry. DiFX has been extensively tested against three existing hardware cor- 
relators and is now used full-time for the LEA National Facility; it is also presently being 
adopted by other VLBI arrays, including the VLBA. The astrometry project has been suc- 
cessfully concluded with the measurement of seven pulsar parallaxes, more than tripling 
the number of Southern Hemisphere pulsars with a VLBI-measured parallax. Thus, the 
twin goals set at the outset of the thesis have been satisfactorily accomplished. The con- 
tributions of this thesis to, and possible future directions of, the main thesis areas of VLBI 
instrumentation and VLBI pulsar astrometry are considered in more detail below. 



7.1 VLBI instrumentation 

The DiFX software correlator developed during this thesis continues to expand in capa- 
bility and scope, and is expected to become the primary correlator for a minimum of two 
additional arrays (the VLBA and the new Australian/NZ geodetic array funded under 
AuScope) during 2009. Unique experiments such as extremely high spectral resolution 
pulsar scintillation studies, using high sensitivity VLBI, have been correlated using DiFX. 
These experiments would not have been feasible using existing correlator infrastructure. 
The first eVLBI observations within Australia, and the first Australian eVLBI observa- 
tions incorporating Asian antennas, have been made using DiFX. As more telescopes are 
connected to high-speed networks, eVLBI offers the possibility of dramatically reducing 
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operating costs for VLBI networks, by reducing costs associated with maintaining and 
shipping large quantities of physical storage. 

For the LBA, the bandwidth upgrade enabled by DiFX has resulted in a significant 
increase in the number of proposals to use the instrument. Other substantial improve- 
ments have been made to LBA functionality during this thesis, such as reductions in the 
uncertainty of antenna positions and refinement of the correlator geometric model. These 
improvements are reflected by the number of ongoing astrometric LBA programs, and the 
number of current astrometric LBA proposals. 

Future developments in DiFX will see the addition of capabilities such as single pass 
correlation of multiple phase centres, adaptation to new recording and output formats and 
increased ease of use. These features, coupled with other instrumental improvements such 
as wide-bandwidth digital backends, will allow existing VLBI arrays such as the LBA 
and VLBA to begin to conduct SKA-style wide-area, high angular resolution surveys 
at high sensitivity. As well as dramatically improving knowledge of the radio sky at 
high resolution, this will help pave the way for SKA observations through development of 
wide-fleld imaging techniques and calibration. 

As the power and efficiency of digital electronics continue to grow in the future, an 
ever-larger region of interferometer parameter space will become economic to correlate 
with software-based, as opposed to hardware-based, solutions. In the three and a half 
years of this thesis, the computing capacity of the Swinburne supercomputer has increased 
four-fold, and most correlations are now run using only the CPUs available on the machines 
which host the storage media. For existing and planned VLBI arrays (excluding the SKA), 
a software-based solution is feasible now, and will become more so in the future. This 
offers the potential to continue progress towards the standardisation of VLBI formats, 
which with attendant improvements to support software and usability will help make 
VLBI observations accessible to a wider cross-section of the astronomical community. 



7.2 VLBI pulsar astrometry 



Of eight pulsars observed in this thesis, parallax and proper motion measurements have 
been made for seven, with a proper motion determination for the final object. This is 
comparable to the largest published N orthern Hemisphere parallax program conducted 
using the VLBA (jBrisken et al.l . 12002), where nine parallaxes were observed from ten 
objects. Ironically, the offending pulsar was the same in both instances - PSR J2145- 
0750. 
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The ensemble of independent distance and velocity measurements for these seven pul- 
sars has been used along with previously published results to analyse the capabilities and 
shortcomings of the NE2001 Galactic electron distribution model used to predict pulsar 
distances, the accuracy of estimates of pulsar transverse velocities using scintillation, and 
the presence and relative contribution of different sources of systematic error in relative 
astrometry. This analysis has shown that while the NE2001 distance model offers signif- 
icant improvements over the previous TC93 distance model, appearing to eliminate the 
systematic bias seen in the latter, it still possesses a long tail of large distance errors. 
A simple model of these infrequent, large distance errors shows they systematically bias 
pulsar velocities derived using the NE2001 model, leading to an 18% increase in the mean 
observed 2D transverse pulsar velocity and a fifty-fold increase in the the number of very 
fast moving pulsars (transverse velocities greater than 1000 km s~^). However, this simple 
model is itself likely to be somewhat biased by selection effects in the available indepen- 
dent pulsar distance measurements, and thus these results are almost certainly an upper 
limit to the errors induced in velocity estimates. An unbiased parallax sample would aid 
in quantifying the true impact of distance model errors. Scintillation velocity estimates 
have been shown to be unreliable for individual pulsars under the typical assumptions of 
a thin scattering screen with an isotropic turbulence distribution. 

In addition, the distance and velocity measurements have allowed a number of impor- 
tant conclusions for the individual objects, with the most significant being: 

• The radio luminosity of PSR J0108-1431, while low, is a factor of 3.4 higher than 
previously thought. Interestingly, though, its recently measured x-ray luminos- 
ity, coupled with its newly measured distance and proper motion implies that it 
is the most efficient converter of spin-down luminosity to x-ray luminosity amongst 
rotation-powered pulsars, and that reheating must occur for isolated neutron stars. 

• The most precise pulsar distance determination made to date, for PSR J0437-4715, 
constrains the time rate of change of Newton's gravitational constant G to be less 
than 3 parts in 10^^, and precludes the presence of an unseen Jupiter-mass planet 
within 226 AU of the Sun (within 50% of the sky). 

• The parallax of PSR J0630-2834 shows that the pulsar is more than four times closer 
than predicted by the NE2001 distance model, revising its x-ray luminosity down 
by a factor of 19 and proving that the system does not require an unusually efficient 
method of producing x-rays. 
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• The first significant measurement of the parallax of the double pulsar J0737-3039A/B 
shows that the system is more than twice as distant as the DAf-derived predictions, 
and enabled calculation of the kinematic terms which contribute to the observed 
orbital period derivative. The accuracy of the distance and velocity determination 
from this thesis will permit tests of GR using the double pulsar system to proceed 
to the 0.01% level. Furthermore, the confirmation of the system's low transverse 
velocity supports formation scenarios in which neither pulsar received a large kick 
at birth. Finally, the revised distance strongly supports a magnetospheric origin for 
the majority of the x-rays observed in the PSR J0737-3039A/B system. 

• PSR J2144-3933 is shown to be even closer than its DM-derived distance of 180 
pc, confirming that it is the least luminous pulsar known at 1400 MHz. The revision 
to its period derivative based on the measurement of distance and proper motion 
places it even further past the traditionally assumed pulsar "death line" . 

The future of VLBI pulsar astrometry is bright, with continual improvement in instru- 
mentation bringing more and more pulsars within reach of the technique. With this work 
having shown that the measurement of pulsar parallaxes with the LEA is possible out 
to and beyond two kpc, a large number of southern pulsars are accessible. Several LBA 
proposals for pulsar astrometry are currently active, and the number of VLBI parallaxes 
in the south should continue to grow. A proposal has also been submitted for the contin- 
uation of astrometry on PSR J0437-4715, which could lead to the most accurate VLBI 
parallax measurement of any object to date, and allow a limit on the time rate of change 
of G which surpasses that presently made by Lunar Laser Ranging. 

Beyond the LBA, the sensitivity upgrade of the VLBA which is presently underway 
should deliver 4 Gbps recording capabilities during 2009. This upgrade will deliver in- 
stantaneous sensitivity surpassing that of the LBA and approaching that of the EVN, but 
with a much greater field of view. When combined with the full-time operation and con- 
current systematic control available with the VLBA, the upgrade should allow the deepest 
VLBI-resolution imaging ever made. This will enable almost any known pulsar to be tar- 
geted for VLBI astrometry, and allow larger and more complete astrometric surveys to be 
undertaken. Further in the future, the SKA will integrate VLBI into connected-element 
observing and make astrometry of pulsars almost routine; it is likely to revolutionise our 
knowledge of their spatial, velocity and luminosity distributions. The trickle of mea- 
sured pulsar parallaxes is likely to become a flood over the next decade, and may help 
astronomers finally gain a solid understanding of the birth, life and death of pulsars. 
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